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ABSTRACT 


This report documents a study which was carried out using a 
Monte Carlo Simulation method to determine the effect of scattering 
from turbid water on the polarization of a backscattered beam of 
laser light. The relationship between the polarization and the type 
and amount of suspended particulates in the water was investigated. 
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l.o INTRODUCTION 


Excessive amounts of suspended particles in the environmental 
waters lead to significant economical and ecological consequences. 
Economic losses are mostly due to expenditures for dredging 
operations and the reduction of the nation's reservoir capabilities. 
Less tangible, but probably more important are the adverse 
ecological effects. For example, organic debris and residues 
absorbed on the surface of particles increase biochemical oxygen 
demand (BOD) and chemical oxygen demand (COD) which cause a concomitant 
decrease in the concentrations of dissolved oxygen and the elimination 
or reduction of sensitive organisms. The introduction of sediments 
into an aquatic environment also may alter nutrient concentrations, 
especially total organic carbon (TOC), nitrogen, and phosphorus. 
Increases in nutrient loading may increase productivity in the 
absence of inhibitory factors. Although high productivity can be 
beneficial, it also may lead to serious eutrophication and nuisance 
algae blooms. As long as the concentration of suspended sediments 
remains high, primary productivity (photosynthesis) will be inhibited 
or much reduced due to a reduction in the penetration of light in 
the water column. Sediments also may transport toxic materials such 
as pesticides and certain heavy metals which may adversely affect 
the biota. In the light of the accumulated evidence, the know- 
ledge of sediment concentration and composition seems of utmost 
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importance to the understanding of ecological processes in the 
aquatic environment. 

Considering the dynamic character of the environmental waters 
monitoring procedures for measuring water quality parameters should 

be based on a timely data collection system, such as can be provided 
by applications of remote sensing technology. To develop appropriate 
remote sensing tools a laboratory field program is presently 
being pursued at the NASA/Langley Research Center (LaRC) . The purpose 
of this program is to investigate the remote sensing of water quality 
parameters using information on the polarization properties of the 
backscattered radiance when a polarized laser beam is directed at 
turbid water. The present report describes a modeling effort which 
deals with variations in the polarization characteristics of the 
backscattered radiance as a function of water turbidity, suspended 
particulate type, and detector spot size. This investigation employs 
Monte Carlo simulation techniques to describe the radiative transfer 
processes in the turbid medium, and to determine the backscattered 
radiance. Specifically the polarization properties of the back- 
scattered beam is calculated based on the multiple scattering events 
and the optical properties of the water medium. 

The organization of the report is as follows: Section 2.0 

presents a summary of the investigation and the major conclusions; 
the details of the METREK Monte Carlo model are given in Section 3.0 
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while the details of the calculation of the scattering functions are 
presented in Section 4.0; Section 5.0 discusses the effect of 
scattering on polarization in an abstract, theoretical framework while 
Section 6.0 discusses the effect of scattering based on the Monte 
Carlo simulation and gives the results of the calculations. 
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2.0 SUMMARY AND CONCLUSIONS 


This work has been done in support of an on-going laboratory and 
field experiment at the NASA/Langley Research Center (LaRC) . The 
goal of this work is to investigate the effect of varying the 
amount and type of suspended particulates in water on the polarization 
of a backscattered beam of laser light. In addition, it is expected 
that the results of the study will provide some guidance to the experi- 
ment design. 

The first step in performing this work is the calculation of 
the optical properties, represented by the Mie matrix, of selected 
soil types. Information concerning the size distributions of these 
soils was obtained from a previous LaRC study . The calculated 
quantities are the Mie matrix and the associated volume scattering 
phase functions and volume scattering distribution functions. It 
is found in this portion of the study that the dominant source of 
variation in the scattering functions is the maximum size of the 
particles included in the size distribution. This has important 
implications for attempts to compare the laboratory results (where 
the water is constantly mixed) and the field test results (where 
the larger particles may -have settled out) . 

A Monte Carlo simulation is next run using one of the calculated 
scattering functions. This simulation is done in such a way so as 
to follow each scattering event of each photon which emerged from 
the water. The results of the 'simulation are then analyzed using 
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a method which allows the calculation of the resultant intensity 
and polarization of the backscattered beam for different detector 
spot sizes (the size of the area on the surface of the water from 
which photons emerging from the water will enter the detector) , 
different concentrations, and different optical properties (in terms 
of the ratio of the amount of absorption to the amount of scattering 
of the suspended sediments). 

Based on the simulation study, the use of a polarization 
measurement of a backscattered beam of laser light for the determina- 
tion of the concentration of suspended particulates appears feasible. 
The usefulness of the technique for discriminating between different 
types of suspended particulates is less certain. In this study a 
change in particulate type was represented by allowing the relative 
amount of scattering and absorption to vary over a limited range. 

It is found that there is a small dependence of the polarization on 
the relative amounts of absorption and scattering over the range of 
interest. In fact a change in particulate type will result in not 
only a change in the relative amounts of absorption and scattering 
but also m a change in the volume scattering distribution function. 
This effect has not been examined in this study. The results of 
this study suggest that a trade off exists between the detector spot 
size and signal- to-noise considerations. The ability to discriminate 
between different particulate concentrations is improved if the 
detector spot size is decreased. Thus the smallest detector spot 
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size commensurate with signal-to-noise requirements is indicated by 
this study. 

While only a limited number of concentrations s optical properties 
and experiment configurations were considered in this study a simple 
phenomenological theory has been developed which should allow an 
extension of the results to other cases of interest. This theory 
relates the polarization of the backscattered beam to the detector 
spot size and the concentration of suspended particulates and gives 
excellent qualitative agreement with the results of the simulation. 

Comparison of the theoretical predictions with data obtained 
from an experiment in the Chesapeake Bay area (see Figure 6-7) shows 
very good agreement. 
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3.0 METREK’S MONTE CARLO RADIATIVE TRANSFER MODEL 


The important optical quantity which represents the scattering 
properties of a medium is the Mie Scattering Matrix, S(m,0), which 
determines the intensity and polarization of the radiation scattered 
into an angle 0 from a volume element containing particles with an 
index of refraction m (see Section 4.0 for a more thorough discussion 
of the Mie Matrix). A number of optical parameters of interest can 
be derived given knowledge of the Mie Matrix. These include the 
absorption coefficient, a, the scattering coefficient, s, and the 
total extinction coefficient, a = a + s. (For a discussion of the 
relation between these parameters and the Mie Matrix see Appendix A) . 
The physical importance of these quantities can be seen from Beer's 
law which states that, for a collimated beam of light passing through 
the medium, the intensity decreases exponentially as a function of 
distance traveled, r: 


-a r -ar 

I ^ e = e 


-sr 


(3-1) 


m. qje 

That is, e measures the fraction of the beam intensity which is 

*”S3T 

lost due to absorption and e measures the fraction of the beam 
intensity which is lost due to scattering out of the collimated beam. 

Two other important optical parameters are the volume scattering 
phase function o(0) and the scattering distribution function F(0) . 
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The scattering phase function o(0) is related to .the Mie Matrix 
through 


o(0) = S 1;L (m, 0) + S 22 (m,0) 


(3-2) 


and represents the total intensity scattered into an angle 0 from 
a volume scattering element ,0properly normalized this is equivalent 
to the probability of scattering into an angle 0). The 
scattering distribution function is the cumulative probability of 
scattering between 0 and 0 degrees and is given by 


f(0) 


J o(0’ ) Sin 0' 

d0 r 

IT 


/ 0 (0 ' ) Sin 0' 

d0* 


(3-3) 


The role of each of these parameters in the Monte Carlo simulation 
will be discussed in the next section. 

3.1 Monte Carlo Simulation for Narrow Beam Transmission 

The Monte Carlo simulation technique is based on the use of a 
random number generator to determine when each photon undergoes a 
scattering event and into which set of angles (6,<j>). For this study 
the photons are assumed to strike the water surface at normal 
incidence in a narrow beam. In addition it is assumed that each 
photon is linearly polarized in the same sense and of the same 
wavelength (A = 500 nm) . Thus the simulation is that of a polarized 
laser beam. 

The photons that are incident on the water surface are 
refracted at the air-water interface in accordance with Snell's law 
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(see Figure 3-1). To decide how far the photon travels before a 
scattering event occurs, a random number r^ is chosen in the uniform 
interval (0,1) and the distance AB is set according to the formula 
AB = - ln(r ) (3-4) 

AB is in units of scattering length, s \ 

At point B two new random numbers, r^, r^, are chosen and used 
to determine the scattering angles 02 , and (©» <j> are measured m 
a sperical coordinate system) . The ij>-angle is chosen to be uniformly 
distributed, i.e., ^ = 2irr^. The selection of 02 is accomplished 
from the relation r 2 = F(0), where F(8) is the scattering distribution 
function for polar angle. This is equivalent to chosing 0 with a 
probability distribution identical to that given by the volume 
scattering phase function, o(0). After the selection of cj^ and 
a random number is chosen for computing the distance traveled, and 
the process will continue until the photon emerges from the water. 
Internal reflection* from the water- air interface is also treated in 
the model. 

At each step In the Monte Carlo process, the angles of scattering 
are chosen with respect to the previous incident photon direction, 
which is in general different from the direction of the original 
z-axis. In order to keep track of the photons with respect to the 
original coordinate system, (and be able to specify its coordinates 
as it leaves the water) coordinate transformations are applied 
after each scattering event. These transformation are: (1) a 
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AI’R-WATER 

INTERFACE 


Z-axis 


FIGURE 3-1 

MONTE CARLO SIMULATION OF NARROW BEAM TRANSMISSION 

IN TURBID WATER 


rotational transformation to set the instantaneous system of the 
photon parallel to the original system, and (2) a translational 
transformation that takes into account the displacement of the 
instantaneous system with respect to the original system. 

The procedure described so far takes into account only the 
scattering processes. In order to take absorption into account 
the probability weight associated with each photon emerging from 
the air-water interface is reduced by a factor of e 1 , where y is 
the total path length traveled in the water by the photon, and a 
is the absorption coefficient. 

For each photon that emerges from the water surface at an angle 

from the normal of less that 10° a record is kept of the 0 and tj> 

angles for each of the scattering events. Thus the statistics 

generated are applicable to the case of having a near nadir looking 

detector. In addition, the point of exit from the water is 
-1 

measured from the entry point in units of s , the inverse scattering 
coefficient. Thus any detector spot size can be accomodated using 
the same set of statistics - all photons whose exit point is less 
than (s x detector’s spot size) can be considered to have reached the 
detector. Clearly this allows the statistics generated from any 
simulation run to be used to analyze any combination of spot size 
and s values. The only parameter which is fixed in the Monte Carlo 
run is the scattering distribution function F(6). The values for 
spot size, scattering coefficient and absorption coefficient, and 
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extinction coefficient, can all be varied after the initial 
simulation. For this reason, values for s, a, and a are chosen 
as adjustable parameters rather than being actuall-y calculated from 
the Mie scattering formalism, by allowing s (or a) to vary while 
keeping the ratio a/s constant we can examine the effect of varying 
the concentration of a particular sediment. By the same token, 
holding s constant and varying the 'a/s ratio allows the investigation 
of the effect of changing the optical characteristics of the sediments. 
For each set of spot size, s value, and a/s value the photons which 
reach the detector are selected from the statistics of the Monte 
Carlo run and their Stokes vectors are computed using the set of 
angles '(0's and <j>'s) that each photon has been scattered. (For 
details see section 5.0). The resulting set of individual photon 
Stokes vectors are then averaged to yield a total Stokes vector for 
the backscattered beam which is received by the detector. Given 
this Stokes vector the total polarization and* the degree of linear 
polarization can be calculated. 
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4.0 CALCULATION OF SCATTERING FUNCTIONS 

This section is devoted to the theoretical treatment of 
scattering and absorption from suspended particulates. The Mie 
theory of light scattering from a single particle is treated in 
Sub-section 4.1. The extension of Mie theory to the case of 
polydisperse suspensions is then discussed along with the computa- 
tional methods used to calculate the scattering function, in Sub- 
sections 4.2 and 4.3 (Appendix A discusses the relationships between 
the Mie parameters and the extinction, scattering and absorption 
coefficients). Sub-section 4.4 includes a discussion of the size 
distributions and optical properties of the clay sediments considered 
in the calculations. Finally, rin Section 4.5 the results of the 
calculation of the scattering functions are presented. 

The following discussion of the Mie theory of scattering and 
the computational methods is a brief summary. For more detailed 
discussions of Mie theory for single scattering, the reader is 
referred to References 2, 3, and 4. Reference 5 contains a 
discussion of Mie scattering from polydispersions and Reference 6 
contains the details of the computational procedures and requirements. 

4.1 Mie Theory for Single Particle Scattering 

When light is incident on a particle, it undergoes both scattering 
and absorption (the inelastic scattering processes x$hich result in a 
change in frequency are ignored in this study) . The characteristics of the 
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scattered radiation depend on the wavelength, X, of the incident 
light, 'the generally complex index of refraction, m, of the particle, 
and the size, r, and shape of the particle. In this report we will 
restrict the discussion to spherical particles; for the treatment of 
inorganic sediments in water this is probably not a serious 
restriction. 

■If a monochromatic beam .of light, represented by the Stokes 

vector I „ is incident on a spherical particle at an angle 0 = 0, 

•o 

‘then the ‘Stokes vector ,o£ the light scattered is given by 


I(x,m, 0) = S(x,m,0) I , (4-1) 

. Z 'v O 

, 4ir - , 

where S'(x,,m, 0.) is the single particle scattering matrix which depends 
in general, on the size parameter 


x = 


2irr 

X 


(4-2) 


and the complex index of .refraction,, m. The calculation of S(x,m, 0) 
requires the solution of Maxwell's equation in spherical coordinates 
with a discontinuous change in the index <of refraction across the 

(7) 

spherical surface. This solution was originally derived by G. Mie 
in 1908, and independently by P. Debye^ in 1909. 
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The scattering matrix can be written as: 


S(x,m, 0) 


/ M 1 (x,m,0) 

0 

0 

° \ 

0 

M 2 (x,m,0) 

0 

0 1 

0 

0 

S 21 (x,m,9) 

-D 21 (x,m,0) i 

0 

0 

D 21 (x,m, 0) 

S 21 (x,m,0)/ 


and the Mie solution is 

M^x.m,©) = S 1 (x,m,0) S 1 *(x,m,0) 

M 2 (x,m,0) = S 2 (x,m, 0) S 2 *(x,m, 0) 

S 21 (x,m,0) = %(S 2 (x,m,0) S-^x.m, 0) + S-^XjitijQ) S 2 *(x,m,0)) 

D 21 (x,m,6) = Js(S 2 (x,m,e) S 1 *(x,m, 0) - S 1 (x,m,0) S 2 *(x,m,0)) 

Where S^(x,m, 6) and S 2 (x,m,8) are the complex amplitudes for the 
scattered radiation (and * designates the complex conjugate). 


S 1 (x,m,0) 


S 2 (x,m, 0) 


Y 1 (2n+l) 

n=l n(n+1> 


(2n+l) 

n ( n+1 ) 



(x,m)ir n (p) + b n (x,m)T n (y) 


(x,m)ir (y) + a (x,m)r (y) 
n n n 


} 


(4-5) 


In these expressions ^(h) and t (yi) are derivatives o;f the Legendre 
Polynomials : 


T n <4) 


T (y) 

n 


dP n (y) 

dy 


PiT n (y) - (1-y > 


dir n (y) 

"dy 


(4-6) 
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(where y = cos 0). Also 


^(nK}ip n (-x-) - nnJJ n (mx)^(x) 

a n (x 3 m) - ^ nn|» n (mx)?^(x) 

(4-7) 

nufi^(mx)^ n (x) - (p n (mx)^(x) 

b n ( 'X’ m) = mtfT (mx) ? n (x) - ^ (rax) S' (x) 

and the ifi's and g’s are related to the spherical Bessel functions of 

the first and second kinds (j^ and respectively): 

il> (z) = zi (z) 

y n ' J n 


? n <x) •= xj n (x)-iy n (x) 

4^(z) = zj n _ 1 (z)-nj n (z) 

? n (x) = X j n _i ( x)-i yil _i(x)-n3 n (x)-iy n ( X ) 


(4-8) 


,4.2 Mie Theory for Scattering from Polydispersions 
A polydispersion is a suspension of scattering particles of 
uniform physical characteristics hut of varying number concentration 
depending on particle size. Because of the existance of different 
particle .sizes it makes little sense to talk of scattering from a 
single particle. Instead, it is useful to consider the scattering 
properties of a small volume element containing a number of particles. 
The size of this volume element is of some, at least theoretical, 
importance. Clearly, if it is to be used to -represent the scattering 
properties of all similar volume elements then it must contain a 
representative set of particle sizes - this requires that the volume 
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element not be too small. On the other hand, since we are consider- 
ing only single scattering from the volume element, it must not be 
too large. An additional condition that must be imposed is that the 
interparticle separation be large compared to the wavelength. The 
reason for this is that the interaction of light with a particle 
will be assumed independent of the interactions with all other 
particles. This condition requires that the particle density in the 
volume element not be too large. For our purposes, it will be assumed 
that all of the above conditions are satisfied. 

The polydispersion can be completely specified, for our purposes, 
by an index of refraction m and a probability density function n(r). 
The density function gives the relative concentration of each size 
contained in a volume element. 

The characteristics of the scattered radiation due to the volume 
element can then be represented by a volume scattering matrix 
£(m, 0) in a manner analogous to Equation (2-1): 

2 

l(m,0) = §(m,0)I o (4-9)' 

4ir 

This volume scattering matrix can be calculated from, the set of 
particle scattering matrices: 

CO 

S(m,0) = J £(x,m, 9) n(r)dr (4-10) 

o 
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where 


CO 

J n(r)dr = N (4-11) 

o 

and N is the total number of particles per unit volume. In what 
folloxtfs, N will be assumed to be unity since S(m, 0) scales with N. 
The ability to represent S(m, 0) as a linear superposition of the 
individual scattering matrices is a direct consequence of our 
assumption that the interparticle separation is much greater than X. 

The calculation of S(m,0) thus reduces to calculation of the 
individual S(x,m,0) and then integration over all sizes with the 
proper weighting given by n(r). 

4. 3 Computational Methods 

The calculation of the scattering phase functions and the 

averaging over size distributions was carried out on an IBM 370/148. 

The program listings are reproduced in Appendix D. 

In computing the sums in Equation (4-5), the major difficulty 

arises in the evaluation of the a n (x,m) and b n (x,m). Using the 

'definitions of ip , \p' a E, , and E,' , and the standard recurrance 
n r n n’ hi’ 

relations for the Bessel functions. Equation (4-7) can be rewritten: 
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where 


A (mx) 
n 


^(mx) 

if> n (mx) 


(4-13) 


the logarithmic derivative of ^(mx), and Re denotes the real part. 
The natural approach to the evaluation of Equation (4-12) is to 
employ a standard (upward) recurrance procedure. Unfortunately, if 
the imaginary part of the index of refraction is not zero and n is 
large then the upward recurrence procedure results in large instabil- 
ities in the calculation of A (mx) . For this reason, the DBMIE sub- 

n 

routine employs a downward recurrence procedure to calculate the 

A (mx)s. These values are then stored for use in the evaluation of 
n 

Equation (4-12). Becuase of the large storage requirements resulting 
from this proceudre (n n. 7000), and the fact that double precision is 
employed in all of the calculations, a virtual machine with 512 K 
bytes of storage is required for the implementation of the DBMIE 
and POLYMIE routines . 

The scattering phase functions are computed in the DBMIE sub- 
routine, and the average, Equation (4-10) , is computed in the calling 


routine POLYMIE. While analytic functions have been used for 'the 
size distributions, n(r), the integral has been approximated by a 
summation over a discrete set of radii. Tests to determine the 


effect of using a summing procedure have shown that this results in 
no loss of accuracy. In addition, test runs were made to compare 
the results when Ar = O.lp and Ar = lp were used in the summing pro- 
cedure. The use of Ar = lp resulted in no significant change in the 
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results from those obtained using Ar = . ly over the range 0 < r < 
lOOg. Calculations were made using r = 100y and r = lOy 

TQ3.X H13X 

(Ar = 0-ly). A discussion of the proper upper limit for r is given 
in Section 4.4.2. 

The amount of virtual CPU time required for these calculations 
is significant and has been a major factor in determining r max and 
Ar. As an example, the calculation of the volume scattering phase 
function for a polydispersion with m = 1.144 - O.Oi, A. = 0.5y, 
r max = anc * = re U u: i- res approximately 26 minutes of virtual 

CPU time. 

4.4 Properties of Clay Samples 

Data on four different clay samples were provided by NASA/LaRC. 
This data consisted of empirical size distribution curves as well as 
brief descriptions of chemical composition. The physical character- 
istics of the clay are discussed in Section 4.4.1 while the size dis- 
tributions are presented in Section 4.4.2. 

4.4.1 Physical Characteristics of Clay Samples 

Four types of clay were selected by NASA/LaRC. These were: 
Feldspar, Calvert, Ball and Jordan. According to the analysis of 
these clays performed by NASA/LaRC the compositions are: 

• Feldspar - Feldspar and Quartz Minerals 

• Calvert and Jordan - Kaolinite and Illite 

• Ball - Montmorilloite, Kaolinite and Illite 
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The real refractive index and chemical components of these 

(9) 

minerals is shown in Table 4.1. For reasons which will be dis- 

cussed in Section 4.4.2, Feldspar and Ball Clay were chosen to be 
included in this study. 

To estimate the index of refraction of the clay samples, we 
take a simple average of the indices of refraction of the components. 
Thus, for both Feldspar and Ball Clay, the real part of the index of 
refraction is estimated as 

ReCm^r) = 1-53 . 

This, of course, is the index of refraction with respect to air and 
we require the index of refraction with respect to water which can 
be obtained by dividing R e ( m ^ r ) by the index of refraction of water 
1.337 (for wavelengths of approximately 500 nm) . 

Thus 


Re(m . ) = 1.144 . 

water 


Estimating the imaginary part of the index of refraction is not 
so straightforward, since direct measurements of Im(m) have not been 
made. Since these minerals have very low conductivity, it is expected 
that the imaginary part of m will be quite small. The imaginary part 
of m has been measured for soil aerosols and has been found to be 


about . 005 (with respect to air) . 


( 10 ) 


For this study two values for 


Im(m) will be used: 

Im(m ) = 

water 


0 , Non-absorbing 

0.005 = 0.004, Weakly-absorbing 

1.337 
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TABLE 4.1 


CHEMICAL COMPOSITION AND INDEX OF REFRACTION 
OF CLAY .CONSTITUENTS 


NAME 

CHEMICAL COMPOSITION 

INDEX OF 

REFRACTION 

Kaolmite 

A1 2 0 ,.2Si0 2 .2H 2 0 


1.56 

Illi'te 

Vl . 5 A1 4 Si 7-6 . 5 Ai l-l . 5°2 (0H) 4 


1.54 

Montmorilloite 

(..5Ca,Nai) ? (Al,Mn,Fe)^(Si, Al) g ) 

20 (HO) 4 nH 2 0 

1.48 

Fel'dspars : 

Mi crocline 

K 2 0.Al 2 0 3 .6Si0 2 


1.52 

Andes ine 

(030^20) Al 2 0 3 .4Si0 2 


1.55 

Antho class 

(Na,K) 2 0.Al 2 0 3 .6Si0 2 


1.53 
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4.4.2 Particle Size Distributions 


Emperical cumulative size distributions for the four samples 
were provided by NASA/LaRC and are shown in Figure 4-1, 4-2, 4-3 
and 4-4. It is apparent from these figures that the size distributions 
for Ball, Jordan, and Calvert differ significantly from the size 
distribution for Feldspar. Since it was planned that two distri- 
butions would be employed, Feldspar and Ball Clay were chosen. 

This choice allows the investigation of the effect of radically 
different size distributions. 

To utilize the size distribution information, it is necessary to 
determine the size distribution density function, n(r), which specifies 
the relative number of particles with radius r per unit volume. If 
we denote the cumulative size distribution as provided by NASA/LaRC 
as N ( r o > then the relationship between N(r Q ) and n(r) is given by: 


N(r o ) = 


1 - 



dr, 


(4-14) 


or 


n(r) 


dN(r ) 

o' 

dr 


r 

o 


r 


(4-15) 


° I 

A general curve fitting routine (See Appendix E) was used to deter- 
mine the best distribution for both the Ball Clay and Feldspar. 
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Coarser by Weight 
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Particle Size, Millimeters 


FIGURE 4-1 

CUMULATIVE SIZE DISTRIBUTION OF FELDSPAR SAMPLE 




FIGURE 4-2 

CUMULATIVE SIZE DISTRIBUTION OF CALVERT SAMPLE 




FIGURE 4-3 

CUMULATIVE SIZE DISTRIBUTION OF BALL SAMPLE 



Percentage of Particles 
Coarser by Weight 





A plot of this size distribution, density function is shown in Figure 
4-5, while a plot of the corresponding cumulative size distribution 
function (as obtained from Equation 4-16) is shown in Figure 4-6. 

As can be seen in Figure 4-6, the modified Gamma distribution gives 
a good fit to the data points obtained in the NASA/LaRC analysis. 

To fit the size distribution of the Ball Clay sample, Junge's 
distribution model was chosen: 


n(r) = a r 

with the parameters, 

a 1 = 0.2006 

a 2 = 1.624746 


(4-17) 


determined using the same curve fitting routine employed for Feldspar. 
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n(r) - Relative Number of Particles per Unit Volume 



Particle Size, Microns 


FIGURE 4-5 

PARTICLE SIZE DENSITY FUNCTION FOR FELDSPAR (MODIFIED GAMMA DISTRIBUTION) 
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FIGURE 4-6 

CUMULATIVE SIZE DISTRIBUTION FIT OF FELDSPAR SAMPLE USING MODIFIED GAMMA DISTRIBUTION 




FIGURE 4-7 

PARTICLE SIZE DENSITY FUNCTION FOR BALL CLAY (JUNGES DISTRIBUTION) 
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The size distribution density function and the cumulative size dis- 
tribution function for Ball Clay using Junge's distribution are shown 
in Figure 4-7 and 4-8. It is apparent from Figure 4-7 that Junge's 
distribution function is not, strictly speaking, a probability dis- 
tribution since the integral (Equation 4-11), 



dr = N 


can not be normalized, i.e., N is infinite. However, Junge's distri- 
bution has been found to accurately represent particle sizes of ocean 
sediments. In addition, the lower and upper limits of integration 
in Equations (4-11) and (4-10) are not set equal to zero and infinity, 
in practice, allowing Equation (4-11) to be normalized. 

The question of the proper upper limit for Equation (4-10) and 

Equation (4-11) is of more than theoretical interest. From the 

empirical size distributions provided by NASA/LaRC, it appears that 

an upper limit in Equation (4-10) should be chosen as 100 microns 

( 12 ) 

(ym). However, as can be seen in Table 4.2 the settling rate 
for 100 pm particles is on the order of thirty seconds. Thus, the 
history of the particulates in the body -of water is important. If 
the particulates have been allowed to settle, then the size distri- 
butions determined before the particles are introduced into the water 
are inappropriate. In order to investigate the effect of settling, 
two upper limits, 100 ym and 10 ym, were chosen for the integrals of 
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TABLE 4-II 


SETTLING VELOCITIES OF SAND AND SILT IN STILL WATER 


( Source • Amer. Water Works Assoc.) 

[Temperature 50°F, all particles assumed to ha\<e a specific gravity of 2.651 


Diameter of 

particle Order of Size 


Settling Time Required to 

Velocity Settle 1 Foot 


mm. 


mm./sec. 


10 0 
i o 
0 8 
06 
0.5 
0.4 
03 
0 2 
0.1 S 
0.10 
0 08 
0 OS 
0 05 
0.04 
0.03 
0 02 
0 015 
0 010 
0 008 
0 006 
0 005 
0 004 
O 003 
0 002 
0 0015 
0 001 
0 0001 
0 00001 



Gravel 


Coarse Sand 


Fine Sand 


Silt 


Bacteria 
Clay Particle* 
Col'oldal Particle* 


1.000 

100 

83 

63 

53 

42 

32 

21 

IS 

8 

6 

38 
2.9 
• 2 1 
1 3 
0 62 
0 35 
0.154 
0 098 
0 065 
0 0385 
0 0247 
0 0138 
0 0062 
0 0035 
0 00154 
0 0000154 
0 000000154 


0 3 seconds 
3 0 seconds 


38.0 seconds 


33 0 minutes 


55 0 hours 
230 0 days 
63 0 years 


36 


ORIGINAL PAGE IS 
OF POOR QUALITY 




Equations (4-10) and (4-11). Equation (4-11) was used to properly 
normalize Equation (4-10) with respect to the choice of upper limit. 

4. 5 Results of Computations 

The results of the computation of the volume scattering phase 
functions (4.5.1) and the volume scattering distribution functions 
(4.5.2), using the size distributions of Section 4.4, are presented 
in this section. In addition to examining the effect of settling 
on the calculations, the wavelength dependence of the scattering 
functions are also investigated. 

4.5.1 Volume Scattering Phase Functions 

The computed volume scattering phase functions are shown in 
Figures 4.9 through 4.14. 

Figures 4.9 and 4.14 display the extremely large forward scat- 
tering peak which is primarily the result of including the large 
(VLOO ym) particulates in the size distributions. Both the Feldspar 
and Ball Clay phase functions show considerable difference between 
the non- absorbing and absorbing cases at large angle. While it is 
not evident in the figures, the forward scattering peak is larger 
for the absorbing case at small but non-zero angles (9 — 0.5°). 

Figures 4-11 and 4-12 demonstrate the effect of cutting the size 
distributions off at 10 ym instead of 100 ym. The relative size of 
the forward peak is reduced and the difference between the absorbing 
and non-absorbing cases at large angles is reduced. It is inter- 
esting to note that, although the shape or the Feldspar and Ball 
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FIGURE 4-12 

VOLUME SCATTERING FUNCTIONS FOR BALL CLAY (IOjuM CUTOFF X = 500NM) 
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FIGURE 4-13 

VOLUME SCATTERING FUNCTIONS FOR FELDSPAR {lO^M CUTOFF A = 600NM) 
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FIGURE 4-14 

VOLUME SCATTERING FUNCTIONS FOR BALL CLAY (lO^M CUTOFF X = 600NM) 






Clay size distributions are very different, the upper limit on the 
size appears to be much more important in terms of the difference in 
phase functions. 

Figure 4-13 and 4-14 show the scattering phase functions computed 
for \ = 600 nm (with a 10 ym cutoff) instead of X = 500 nm as in 
Figures 4-11 and 4-12. It can be seen that the phase functions are 
not heavily wavelength dependent. In fact, it can be shown that for 
a uniform size distribution and upper and lower limits of zero and 
infinity in Equation (4-10), the volume scattering phase function will 
be strictly independent of wavelength. 

4.5.2. Volume Scattering Distribution Functions 

While the volume scattering phase function describes the angular 
dependence of scattered radiation, a more important function for use 
with the Monte Carlo simulation is the volume scattering distri- 
bution function, F(0), defined by Equation (3-3). The distribution 
function gives the normalized cumulative probability that a photon 
is scattered in the range 0 to 8 degrees. The volume scattering 
distribution functions for the cases considered in Section 4.5 are 
shown in Figures 4-15 through 4-20. 

It is again apparent in Figures 4-15 and 4-16 that there is a 
considerable difference between the absorbing and non-absorbing case. 
The difference due to the Feldspar and Ball Clay size distributions 
is small. 


44 



Probability 
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VOLUM 





As with the scattering phase functions, the use of a 10 ym cut- 
off decreases the difference between the absorbing and non-absorbing 
cases. In addition, the volume scattering distribution functions 
are changed considerably when the 10 ym cutoff is imposed. 

Figure 4-19 and 4-20 demonstrate the small change in the 
volume scattering distribution functions when the wavelength is 
changed. 
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FIGURE 4-19 

WAVELENGTH DEPENDENCE OF VOLUME SCATTERING DISTRIBUTION 
FUNCTIONS FOR FELDSPAR {lO^M CUTOFF) 
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FIGURE 4-20 

WAVELENGTH DEPENDENCE OF VOLUME SCATTERING FUNCTIONS 
FOR BALL CLAY (IOjuM CUTOFF) 



5.0 THE EFFECT OF SCATTERING ON POLARIZATION 


The description of a light wave in terms of its polarization 
properties requires the use of the Stokes formalism. In this 
section the Stokes formalism is developed and the effects of 
scattering on the polarization of a monocromatic light wave is 
investigated. The reader is referred to Reference 13 for a complete 
discussion of the Stokes formalism and polarization. 

5 . 1 The Stokes Formalism 

We consider a monochromatic, polarized electromagnetic wave 
propagating in the z-direction. The electric vector can then be 
written as 


E = A. Cos(wt) 

X JL 

E = A„ Cos(wt+<5) 

y “ 

E = 0 

z 

This wave can be represented using the Stokes formalism: 



(5-1) 


(5-2) 


Where the Stoke parameters are given by 


I 

I 


1 

2 


A 

A 


2 

1 

2 

2 
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.(5-3-) 


,u <= 

2A l A 2 

'Cos 6 

V - 

2A 1 A 2 

Sind 


There are ,a number polarization '"types" possible for .a completely 
polarized beam of light; some examples are given in Table 5-1. 

5.2 The Effect .of 'Scattering on 'Polarization 
If Phe -light beam is not completely polarized then it can be 
broken into two components: a completely polarized component, 

I , and an unpolarized,, or natural component,, I . In terms of the 

p ’ , n 

Stokes vectors of a 'beam of light *these components are 


I = 
-P 


+ jg((i r i 2 ) 2 + u 2 + v 2 )^ 


Js(I 2 -X L ) 4^%(Cl r I 2 )' 2 + U 2 + V 2 ) H 


u 

V 

\ 

- /jsCij+ig) + ^((i r l 2 ) 2 •+ u 2 + vV*! 


hO- x +i 2 ) + hia ^) 2 + u 2 + v 2 ) % 


I N 


(5-4) 


Notetthat if the beam is completely polarized then 


(I 1 +I 2 ) 


(If I 2 ) 


+ U + V 


\ ' 


rn 


* -if* . k"* 
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TABLE 5.X 

STOKES REPRESENTATION FOR SOME POLARIZATION TYPES 


Type 


X 2 

U 

V 

6 


Description of Polarization 

1) 

1 

0 

0 

0 

0 


Linear, Vertical 

2} 

0 

1 

0 

0 

0 


Linear, Horizontal 

3) 

2$ 


1 

0 

0 


Linear, 1st and 3rd Quadrants 

4) 

% 

4 

-1 

0 

0 


Linear, 2nd and 4th Quadrants 

5) 

% 

h 

0 

1 

' ir/2 


Circular, Right Handed 

6) 

h 

h 

0 

-1 

-ir/2 


Circular, Left Handed 

7) 

5/8 

3/8 

3/4 

3/2 

Sin” 1 

4/5 

Elliptical, Right Handed, 1st and 3rd 
Quadrants 

8) 

3/8 

5/8 

- 3/4 

- 3/2 

-Sin” 1 

4/5 

Elliptical, Left Handed, 2nd and 4th 


Quadrants 



(as can be seen from Equation (5-3)) Equations (5-4) reduce to 


I 

P 



I 


N 



(5-6) 


The degree of polarization, P, is defined as the ratio of the intensity 
of the polarized component to the sum of the' intensities of the 
polarized and unpolarized components, 


P 


I P + I N 


(5-7) 


which in terms of the Stokes parameters is 


P = 


[(I-L-Ig ,) 2 + u ? + V 2 


I l + I 2 


(5-8) 


For completely polarized light P=1 while for a light beam which 
containes a natural component P < 1. 

Equation (5-8) can be inverted to give an alternate expression 
for determining if a wave is fully or partially polarized': 


4I lh 1 

[“}: 

f'fully polarized 

L 'I 

u 2 +v 2 

L>J 

Lpartially polarized. 
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This expression will be useful when examining the effect of scattering 
from a spherical particle, or a collection of spherical particles. 

If an incident wave, represented by the Stokes vector 



is scattered by a spherical particle whose scattering matrix is 
given by Equation (4-3) the scattered wave is given by (apart from 
a constant) 

I 

s 

Upon forming ir s for the scattered wave and simplifying it is found 
that 

4l g 1 I g 2 M L (x,m,0) M 2 (x,m,0) 

U 2 + V 2 S 2 (x,m, 0) + D 2 (x,m,0) 

s s /i 

(5-11) 

M^XjH^e) M 2 (x,m, e) 

__ 2 17 • 

S 21 (x,m,0) + D 21 (x,m,0) 


K ih ) 
vW) ' 
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If the wave incident on a spherical particle is completely polarized 
then ir=l and using the definitions of Equation (4-4) one finds 



S 1 ^x,m,0) 

2 

S 2 (x,m, 0) 

2 

IT = 

S 

S 1 (x,m,0) | 

2 

S 2 (x,m, 0) 

2 


= 1 


(5-12) 


and the scattered wave is also completely polarized. This does not 
mean that the polarization type has not been changed. It can be seen 
from Equation (5-10) that the parameters of the Stokes vector are 
altered by the scattering matrix. For example, in general, linearly 
polarized light will be scattered as eliiptically polarized light. 

If, instead of scattering from a single particle, the light 
is scattered off of a volume element containing a collection of 
particles, then the parameters of the scattering matrix must be 
replaced by their averages (denoted by <...>) over the volume 
element. Thus the polarization condition (Equation 5-11) becomes 
^M^(x, m, 0^ ^M 2 (x,m,0}) 


ir = 
s 


^S 21 (x,m,0)} 2 + (n 21 (x,m,0)) 


ir. 


(5-13) 


-Again assuming that the incoming wave is completely polarized (ir=l) 
and expanding the scattering parameters in terms of the scattering 
amplitudes as was done for Equation (5-12) we find 

2 / \ 2 


TT = 

s 


^(xjm,©)) ^S 2 (x,m,e)^ 


^S 1 (x,m, 0) S 2 (x,m,0) X S 1 (x,m,-0) S 2 (x,m,0)) 


(5-14) 
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and using the Cauchy-Schwarz inequality results in 


2 ^S 2 (x,m,e)^ 2 

{s^x.n^e)^) 2 (s 2 (x,ra,0)y 2 


( 5 - 15 ) 


or ir 1. The equality holds if all of the particles in the 
s 

volume element are identical. For a polydispersion the scattered 
light will be depolarized, that is the degree of polarization will 
be reduced. 


5 . 3 Results of Calculations of Polarization for Single Scattering 
When computing the polarization of a beam of light which has 
been scattered from a particle (or collection of particles) it is 
necessary to first select a specific Stokes vector representation 
of the incoming light. For any type of polarization (i.e., linear, 
elliptical) there are an infinite number of Stokes vector represen- 
tations. The representation which is appropriate for a particular 
scattering event depends on the plane of observation (PO) of that 
event. The PO is defined by the direction of the incoming wave 


and the direction of the observed scattered wave. The Stokes vector 


of the incoming wave is written in the representation for which E 
lies in the PO and E^ is perpindicular to the PO. For the purposes 
of calculating the depolarizing effects of single scattering from 
the clay sediments we will use an incoming wave represented by 
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(5-1) 



This is a linearly polarized wave written in a representation which 
does not give a preference to the E^ and E_^ components. In terms of 
the PO this Stokes vector represents a T^ave whose electric vector is 
at a 45° angle to the x and y coordinates. 

The resulting polarization from the scattering of the wave from 
the Eeldspar and Ball Clay polydispersions is shown in Figure 5-1 and 
5-2. While the angular polarization properties of the two polydisper- 
sions differ there is a much more marked difference between the 
polarizations for the non-absorbing and weakly absorbing' cases. 

When a 10 pm cutoff is imposed on the size distributions 
the differences in the polarization for the non-? absorbing and weakly 
absorbing cases are reduced (Figures 5-3 and 5-4). In addition the 
differences between the polarization curves of the Feldspar and 
Ball Clay samples are reduced. 

A change in wavelength, from 1=500 NM to 1=600 NM, results in 
only small changes in the polarization curves. (Compare Figures 
5-3 and 5-5 and Figures 5-4 and 5-6). 

5.4 Depolarization Due to Multiple Scattering 

It is not possible to compute the depolarization due to multiple 
scattering without actually performing a Monte Carlo simulation. 
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FIGURE 5-1 

POLARIZATION FOR FELDSPAR (A= 500 NM) 
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POLARIZATION 
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FIGURE 5-2 

POLARIZATION FOR BALL CLAY (A = 500 NM) 
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FIGURE 5-3 

POLARIZATION FOR FELDSPAR (lO^M CUTOFF X = 500 NM) 
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POLARIZATION 



FIGURE 5-4 

POLARIZATION FOR BALL CLAY (IOjuM CUTOFF \= 500 NM) 


6 ‘ 



POLARIZATION 
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FIGURE 5-5 

POLARIZATION FOR FELDSPAR (10piM CUTOFF X = 600 NM) 
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POLARIZATION 



FIGURE 5-6 

POLARIZATION FOR BALL CLAY (10 pM CUTOFF A = 600 NM) 
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One of the reasons for this is that the polarization of a wave 
represented by a Stokes vector, depends on the plane of observation 
(PO) which for multiple scattering is defined in terms of the incoming 
direction and the direction to the next scattering particle (or 
volume element). Thus, in general, for each scattering event the 
Stokes vector will have to be rotated into the representation of 
the new PO. For a clockwise rotation of the PO about the propogation 
direction, through an angle <f>, the Stokes vector will have to be 
rotated by operating on it with a rotation matrix R: 


i' = 


(5-16) 


Where 


«<♦> - 


For the multiple scattering of a beam the product 


/ Cos^tfi 

Sin% 

%Sin2$ 

°\ 

2 

Sin <j) 

2 

Cos <j> 

-%Sin2(fi 

0 

-Sin2<j> 

Sin2ij> 

Cos2<j) 

0 

1 ° 

0 

0 

l / 


(5-17) 


n 


ip = n s(e.) R(<j>.)i 

1=1 


(5-18) 


must be calculated. In Equation (5-18) 1^ is the final Stokes vector, 
I q is the initial Stokes vector, £(§..) is the Mie matrix for angle 
and the product is over all scattering events, i=l, n, that results 
in the beam being scattered into the detector. 
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In the next section we examine the results of using Monte 
Carlo simulation techniques to investigate the polarization effects 
due to multiple scattering 
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6.0 RESULTS OF MONTE CARLO SIMULATION 


The Monte Carlo simulation was run for two cases, i.e., for two 
different volume scattering functions, F(6). The cases chosen were 
for Ball Clay with a 10pm cutoff at a wavelength of 500 nm and with 
indices of refraction M = 1.144 - O.Oi and M = 1.144 - 0.0041. 

4 

For the former (Im(m) = 0) the simulation was run in 10 photon 

4 

increments until, after 3 x 10 photons had been run, 335 photons 
were found to have emerged from the water. When a small imaginary 
component is added to the index of refraction the forward scattering 
was increased and it was found that after 10 photons only 35 had 
emerged from the water. It was decided that too large a sample of 
photons would have to be run to obtain a significant number of photons 
and only the M = 1.144 - O.Oi case was investigated in detail. 

Before discussing the results of the simulation a phenomenological 
theory is develope d whi ch is useful in interpreting the simulation 
results. 

6. 1 Phenomenological Theory of Depolarization 

As was discussed in Section 5, each "scattering event from a 
volume element of a polydispersion results, in general, in a decrease 
in the degree of polarization. In addition, each scattering event 
results in a change in the type of polarization. Thus the more times 
a photon is scattered the more its degree of polarization is decreased. 
By the same token, if a photon starts out with a linear polarization, 
the more scattering events it undergoes the more it "loses its 
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memory" of its polarization state. Thus, after a large number of 
scattering events, the degree of polarization should approach zero, 
and the ratio of the y- component of the electric vector to the 
x-component, R = Ey/E z , should approach unity. 

To develop an understanding of the relationship between 

depolarization and the parameters - spot size, s, and a values - a 

relationship between them and the average number of scattering events 

must be investigated. Assuming, for the present, 'that a/s is constant 

and that s is a variable parameter (i.e., the concentration of the 

sediments is variable) this relationship is shown in Figure 6-1. 

In a) the case of high concentration (large s) and small spot size 

is represented while b) illustrates the case for low concentration 

(small s) and' large spot size. These two cases are equivalent 

because they will have the same average number of scattering events. 

In fact, the product (spot size x s) is directly correlated with the 

average number of scattering events. Figure 6-2 shows the results 

of the Monte Carlo simulations in terms of the average number of 

scattering events vs. (spot size x scattering coefficient). (Note 

that, for the M = 1.144 - O.Oi case, the points lie along a smooth 

line and in addition the points are virtually identical for the 10^ 

4 

photon and 3 x 10 photon runs). For the M = 1.144 - 0.004i case the 
points show significant scatter (thus indicating the lack of 
statistical significance for this run). Figure 6-1, c indicates that 
by decreasing the concentration (s value) while keeping the spot size 
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FIGURE 6-1 

RELATIONSHIP BETWEEN SPOT SIZE, DENSITY OF SEDIMENTS 
AND NUMBER OF SCATTERING EVENTS 
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Mean # of Scatterings 



Spot Size x Scattering Coefficient 


FIGURE 6-2 

MEAN NUMBER OF SCATTERING EVENTS FOR BACKSCATTERED 
PHOTONS VS. SPOT SIZE x SCATTERING COEFFICIENT 



constant results in some of the photons which undergo a large number 

of scattering events being out of range of the detector when they 

emerge from the water surface. By the same token increasing the 

concentration (s value) while holding the spot size constant, 

results in a photon which undergoes a larger number of scattering 

events still being in the range of the detector when it emerges 

from the water (Figure 6-1, d). 

Based on our previous discussion of the effect of increased 

scattering on depolarization we can now state the hypothesis : 

As the s value (concentration) is increased (for a 
given spot size) the degree of polarization of the 
backs cattered beam will decrease (the R-f actor will 
approach unity) ; as the spot size is increased (for a 
given s value) the degree of polarization will 
decrease (the R-factor will increase towards unity). 

If s is held constant and a/s is allowed to vary then those 

photons which travel a longer distance, and hence have been scattered 

more, will be weighted less as a/s is increased. Thus one would 

expect the degree of polarization would increase as the a/s ratio 

increases. 

A word of caution concerning the phenomenological theory of 
this section is in order. An important factor has been ignored in 
this development, viz., the variation of the depolarization, due to 
single scattering, with the scattering angle. However, based on 
the fact that most of the scattering is in the forward direction, 
where <fhere is little difference in Figures 5-1 through 5-6 it may 
be reasonable to expect any differences at angles greater than 
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40 degrees to have little impact. As will be seen in the next 
section, the results of the depolarization calculations seem to fit 
the hypothesis concerning the relationship between depolarization, 
spot size, s-value and a/s-ratio. 

6.2 Depolarization Calculation From Monte Carlo Results 

The results of the degree of polarization calculations for the 
backs cattered beam from the Ball Clay, M = 1.144 - O.Oi case are 
shown in Figure 6-3. It is found that the degree of polarization 
does decrease assymptotically as the scattering coefficient is 
increased. It is also found that the degree of polarization is 
approximately the same for combinations of spot size and scattering 
coefficient which yield the same product. Another result of the 
simulation is that small spot size will provide better discrimination 
at bigh levels of scattering coefficient. For discrimination between 
different a/s ratios it is found that a large spot size is necessary. 
From Figure 6-4 it can be seen that there is virtually no dependence 
of the degree of polarization on the a/ s ratio for the 1 inch 
radius spot size. 

The results of the calculations of the R-factor are shown in 
Figure 6-5. Here again it is found that the results agree with the 
hypothesis of Section 6.1. That is, for all spot sizes the R-factor 
increases for large s- values. In addition it is found that the 
R-factor has approximately the same value for corresponding pairs 
of spot size and s-value. As with the degree of polarization it is 
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found that small spot size gives the best results for discrimination 

at high values of s. The R-Factor is also found to be relatively 

insensitive to variations in the a/s ratio (from .1 to .2) for 

constant s-value and small spot size (Figure 6-6). 

Recently data on the depolarization of a backs cattered beam 

of laser light has been collected from a portion of the Chesapeake 
( 16 ') 

Bay v 1 . The results of that experiment, which used a 2.5 inch 
radius spot size, is compared with the theoretical curves in 
Figure 6-7. The agreement is obviously quite good. 
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FIGURE 6-7 

COMPARISON OF THEORETICAL AND EXPERIMENTAL VALUES FOR THE R-FACTOR. 
DIAMONDS REPRESENT DATA COLLECTED FROM THE CHESAPEAKE BAY USING 

A Th INCH RADIUS SPOT SIZE. 




APPENDIX A 


RELATIONSHIP BETWEEN EXTINCTION, SCATTERING, AND 
ABSORPTION COEFFICIENTS AND THE MIE PARAMETERS 
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The extinction (a), scattering (s), and absorption (a) coeffi- 
cients for suspended particulates can be calculated using the Mie 
formalism. Using the Mie parameters, a^(x,ra),and b n (x,m) of 
equations (3-7) the extinction coefficient is given by: 


2ir 


/ " 

(2n+l) -j^Re(a n (x,m) + Re(b n (x,m) ) J. n 


(r)dr 


where n(r) is the particle size distribution function and x = 2iTr/A. 
The expression for the scattering coefficient is: 


2 r co 

-k J n?i (2 " +1 > {K 


(x,m) 


+ 


b (x,m) 
n 




(r)dr 


The absorption coefficient is the difference between a and s, thus 

,2 


a = 


2ir 


/ £ <2n+1) {^n 


(x,m)) + Re(b n (x,m)) - 


a (x,m) 
n 


b (x,m) 
n 




(r)dr. 


The values for a, s, and a used in the Monte Carlo routine were not 
calculated in this way because the values explicitly depend on the 
concentration through n(r). Instead a, s, and a were chosen to 
correspond to physically observed values. 

The absorption coefficient depends on the imaginary part of 
the index of refraction, but in a non-trivial way. If Im(m) = 0 


then it can be shown 


(25) 


that 
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j a n (x,m) - h 2 = % 
b n (x,m) - h 2 = h 

Expanding equation (A-4) leads to 

[ Re ( a n (x,m)> ] 2 - Re(a n (x,m)) + [lm(a n (x,m) ) ] 2 + h> = h 
or 

Re(a n (x,m)) = [Re(a n (x,m))] 2 + [lm(a n (x,m) ) ] 2 
= a n (x,m) 2 

with a similar result holding for b n (x,m). Using these results in 
equation (A-3) leads to a=o. Thus if the imaginary part of the 
index of refraction is zero the absorption coefficient is also zero. 

If Im(m) ± 0 then 

I 2 

a n (x,m) - h < h 

i 2 

b^(x,m) - % < h 

Which, after expansion, leads to 

2 

Re(a n (x,m)) > a (x,m) 

2 

Re(b n (x,m)) > b (x,m) , 

so that, by equation (A-3), a>o for a non-zero imaginary component 
in the index of refraction. 


83 



appendix b 


listings for MONTE’ 


CARLO ROUTINE 


Preceding hage wot ffefvi&u 
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IV G1 RELEASE 2.0 


MAIN 


DATE a 77264 


13/25/51 


C 

c 

c 

C MONTECARLO PROGRAM WITH DOCUMENTATION 

DIMENSION FIS 1100), TESI 100), VALU (501. ANGU50J 
READ 15,25 ) MAXNPH, N MAX, I S 

25 FORMAT 13C2X,I8)I 
READI5.30) TETAI,FI1 

30 F0RMATI2C5X,F8.3)J 
HRITE (6,26) MAXNPH 

26 FORMAT 1*0* , ‘MAXIMUM NO. OF PHOTONS TRACED=‘,I8) 

WRITE (6,27) NMAX 

27 FORMAT ( *0 ', 'MAXIMUM NO- OF EVENTS FOR EACH PH0T0N=‘»I8) 

WRITE (6,29) IS 

29 FORMAT ( '0 INITIAL SEEO FOR RANDOM NO GENERATOR^' , 18) 

WRITE (6,31) TETAI.FII 

31 F0RMAT{»0*, ‘INITIAL IETA IN 0=‘ ,F8. 3, • INITIAL FI IN OEG=‘, 

1F8.3) 

C MAXNPH IS THE MAXIMUM NO OF PHOTONS TRACED 
C NMAX IS THE MAXIMUM NO OF EVENTS FOR EACH PHOTON 
PI=3. 141592654 
RNW= 1.334 * 

C RNW IS THE REFRACTION INDEX OF MATER 

DTRC=PI/ 180.0 

C DTRC IS DEGREES TO RADIANS CONVERSION. 

TETAI=TETAI*OTRC 
FI I = F 1 1 * DTRC 
C START THE CALCULATION 

C NPH IS DEFINED AS THE NO OF PHOTONS AT A GIVEN TIME 

NPH^O 

REA0(9»5004) ( ANGL(I),VALU( I), 1=1, 41) 

5004 FORMAT! F 10. 4, El 5.6) 

C RECORDS NO OF PHOTONS TRACED 

10 IF 1NPH.EQ.HAXNPH 1GC TO 2000 
C TEST FOR END OF COMPUTATIONS 
C PHOTON ENTERS THE MEDIUM AT X=Y=Z=0 

C PHOTON ENTERS AT ANGLS TETAI, FII 

TETA=TETAI 
FI=FII 
X=0. 

Y=0. 

2 = 0.000001 

C DECIOE HOW FAR PHOTON TRAVELS BEFORE AN EVENT OCCURS 
RHOO=RANDNQ( IS) 

T=-ALGG(RHOO| 

GAMA=T 

C T IS THE DISTANCE IN 4-1 UNITS PHOTON TRAVELS TO THE EVENT 
C PHOTON IS AT 
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X = X+T*SIN!TETA)*CGS(FI) 

Y = Y+T*SIN!TETA)*SIN(FI ) 

Z = Z+T$COS (TET A) 

GG TO 150 
100 NPH=NPH+l 

C A8S0RPT ION HAS QCCURED.OR PHOTON FAS CCME GUT OF WATER 

C START A NEW PHOTON 

GG TC 10 
150 CONTINUE 

IF (Z) 400,500,500 
400 XINT=X-Z*TAN(TETA)*COS(FI) 

YINT=Y-Z*TAN!T£TA)*SIN!FI) 

WRITE!6,1083 J 

108 FORMAT!* *, • J= • , 18) 

DINT-SQRI<XINT#*2+YINT**2) 

IF !RNW*SIN!TETA).GT.1.0) GO TO ICO 
TETAAR=A&SIN(RNW*SINITETA) ) 

IFITHTAAR -GT. 0.4IGC TO 100 
WRITE {6,4101 0 INT , T ETAAR 

410 FORMAT!* * * * DISTANCE FROM AXIS=* »F8.5,5X, 'POLAR ANGLE=* ,F8.5) 
WRITE (6,4201 FI,XINT,YIN7 

420 FORMAT (' *,' AZIMUTH ANGLE= ' ,F8. 5 , 5X, • XINT= • ,F8. 5,5X, * YINT=* ,F8.5 
CTA=COS tTETA) 

ACT=A8S ICTA ) 

TCUT=IABS(ZRI— ABS(Z))/ACT 
GAM A=GAMA+T CUT 
WR ITE ( 6 , 109) GAMA 
WRITEI4, 60021 GAMA 
WRlTfc(8,6002) GAMA 
WR l TE !4 , 6002) DINT 
WR ITE 1 8, 6002 )' DINT 
6002 FORMAT I E12.61 

109 FORMAT ( *0' , 'GAMA = *,F8.5) 

WRITEI6, 101) NPH 

101 FORMAT! *0', 'NO. OF PHOTONS TRACED = *,'18) 

JJ=J-1 

WRITE!4, 6001 ) JJ 
WRITE(8,6001) JJ 
6001 FORMAT (15) 

DG 5001 111=2, J 

TES(III)=TES(III)/DTRC 

WRI TE (4, 5003 ) FIS! II I) ,TES(III) 

WR ITE! 8 , 5003 ) F IS { 1 1 1) ,TES ! 1 1 1 ) 

5003 FORMAT { 2F12.6) 

" 5001 WRITEI6, 5002)111, FIS (1 1 1) , TES C 1 1 1 ) 

5002 FORMAT (5X, 13, 5X, • FI= * , F8-3, 5X , ' TETA= • , F8.3) 

GC TO 100 
500 CONTINUE 
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I IV 01 KCLtASE 2.0 


MAIN 


CATE = 77264 


13/25751 


C 2>0 TH15 MEANS PHOTON STILL IN WATER 

00 1250 J=2 * NMAX 

WHAT ARE THE COORDINATE OF TFfc ENC POINT IN THE NGN-ROTATED SYSTE 
FIRST STEP IS TO ROTATE THE SYSTEM USING TETA AND FI 
ROTATION MATRIX IS AIJ, FOR I=J=l-3. 

GENERATION OF THE RLTATIQN MATRIX, WITH THE CONSTRAINT THAT Y- 
AXIS LIES IN A PLANE PARALLEL TO THE YZ-PLANE 
CT=CC S ( ThT A ) 

CF=CC S ( F I) 

CT2=CT*CT 
CF2=CF*CF 
ST= SI N ( T ET A) 

SF=S1MFI> 

ST2=ST*ST 
SF2=SF*SF 
SSL=CT2+SF2*ST2 
SS=SQRT ( SSI J 
SSO=I./SS 

AU=SCRTU.-CF2*ST2 3 
AL2=-Sf-*CF*ST2*SSD 
A13=-CT*ST*CF*5SD 
A22=CT*SSD 
A23=— SF^ST^SSD 
A31=CF*ST 
A33=CT 
A 32=S P*-S7 

ROTATION MATRIX HAS BEEN GENERATED 
SCATTERING HAS OCCURED 

CALL ANGELS FIP, TETAP TO DISTINGUISH FRGM FI, TETA 
FIP* TFTAP ARE DETERMINED IN SYSTEM WITH 2-AXIS 
PARALLEL TO THE INCIDENT DIRECTION 
RHOF=RANDNO (IS) 

FIP=2.*PI*RH0F 
RHGT=RANDNO I IS) 

C PROBABILITY SCATTERING FUNNCT ION FCR POLAR ANGELS FOLLOWS 
DC 1001 1=1,41 

IF (REOT .GE. VALUII3 .AND. RHOT -LE. VALUI I+I) ) GOTO 1002 

1001 CONTINUE 
GO TO 2000 

1002 TETA=ANGLI I ) + ( ANGH I +1)-ANGL I I ) ) * ( RHOT-VALUI I ) ) / l VALUl I+D-VALUI 
ID) 

C CONVERT TETA TO RADIANE 

1011 TCTA=7ETA*DTRC 
TETAP=TETA 

C HOW FAR BEFORE AN EVENT OCCURS, IN THE ROTATED SYSTEM 

FISI J1=FIP 
TESUJ = TCTAP 
RHCJD=RANDNO( IS) 
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T=- ALGG ( PROD 1 

C CALCULATE COOPUINATE OF END POINTS IN ROTATED SYSTEM 

XSTAR=T*SIN[TETAP)*CCSEFIP) 

YSTAR=T*SINITETAPJ*SIN[F1P) 

Z5TAP=T*CQSITETAP) 

C APPLY ROTATION MATRIX TO DETERMINE THE COORDINATE OF THE 

C END POINT IN A SYSTEM PARALLEL TO THE ORIGINAL ONE BUT 

C DISPLACED 

XR=AU*XSTAR+A3l*ZSTAR 
YR=A12*XSTAR+A22*YSTAR+A32*ZSTAR 
2R=A13*XSTAR«-A23*YSTAR + A33#ZSTAR 
C CALCULATE TETA , AND FI IN THE PRESENT SYSTEM, WHICH IS 

C PARALLEL TO THE ORICINAL ONE. 

FI=A7AN1 ABS1YRJ/ABSIXR) ) 

IF (XR.LT.O.O) GO TC 133 
EFtYR) 333,333,6 33 
333 FI=2.*PI-FI 

GC TO 533 
633 F I=F I 

GO TO 533 

133 IF(YR) 233 i 233,433 

233 FI = F I + P I 

GC TC 533 
433 F I = P 1— F I 

533 CONTINUE 

XK2=XR<=XR 
YR2=YR*YR 
ZR2=ZR*ZP 
DT=XR2+YR2+ZK2 
SQDT=SGRT ( OT ) 

TETA=ARCOS < ZR/SGDT) 

C CALCULATE X,Y,Z OF THE END PGINT OF THE PHOTON IN RESPECT TO THE 

C ORIGINAL AXES 

X=X+XR 
Y=Y+YR 
Z=Z+ZR 

IF IZ) 400,400,700 
700 GAMA=GAMA+T 

1290 CONTINUE 
GO TO TOO 

20C0 WRITE (6,5000) IS 

5000 FORMAT {■ * , ’ LAST RANDNO USED=',I10) 

CALL EXIT 
END 
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MAIN 


DATE = 77264 


13/25/51 


C 

C 

FUNCTION RANDNCI IX ) 

I Y= I X*65 535 
IFIIY) 5,6,6 

5 IY=IY+2X47483647+i 

6 KANDNG=IY 

R ANDNC=RANDNO*. 46566 13 E-9 

IX=IY 

RETURN 

END 
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APPENDIX C 


LISTINGS FOR STOKES ROUTINE USED TO 
CALCULATE THE STOKES PARAMETERS OF THE 
BACKSCATTE RED RADIANCE 
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N IV 61 RELEASE 2.0 


MAIN 


DATE 


77270 


13/59/45 


ICO FORMAT ( 15} 

101 FORMAT I 4F5. 2 ) 

104 FORMAT I 2F12.6) 

110 FORMAT! //T10, ' STOKES VECTOR FOR PFGTON NO. ',15} 

111 FORMAT ( 4E12.6) 

113 FCRMAT1//T10, 'FINAL STCKES VECTOR') 

114 FORMAT ( 4E12 .6) 

95 FORMAT! /UO, 'POLARIZATION = ',F7.4) 

REAL*8 TFETA(4.4),SF1NT(500»4) ,FINT(4) , ROT! 4, 4) .DINT 14) .PHI 1500) , 
1TETA ! 20 } ,SMAT!50,6) , STET A (42 ) , SSF IN ! 4 } ,SSSF I! 4 } , GAMA , S ( 20 ) ,SRA! 20) 
1 . RAD (500) 

READ14.150) AS 

READ 14.100) NS 

READ [4, 100) NRAD 

RE AD ( 4. 150 I (SID. 1= l.NS ) 

READ (4 . 150 ) (SR AID . 1 = 1. NRAD) 

150 FORMAT ( 012.6) 

DC 120 J=1 . 4 1 

READ (8, 121) STETA(J), ( SMAT C J . K ) , K=1 , 6) 

120 CONTINUE 

121 FGRMAT( F10.4.5D12.6, F10.4) 

READ (4. 101 ) ( DINT ( I ) *1=1 »4) 

WRITE(6, 130) 

130 FCRMAT(//T10, ' INITIAL STOKES VECTOR') 

WRITE(6,101) (0 INT ( I ) » 1=1 ,4) 

WRITE16.131) AS 

131 F0RMATI//T10. « A/$ VALUE =»,F6.3) 

OG 124 1=1.4 

DC 125 J=l»4 
THETAd, J) = O.ODO 
125 CONTINUE 

FINTI I)=0.000 
124 CONTINUE 
L=1 

102 kEAD(5. 150. END=501) GAMA 
READ 15.150) RAD ( L ) 

KEAD15.100) NSC 

READ (5, 104 J (PHI! J ) , TETA ( J ) , J= 1 , NSC) 

DC 105 1=1,4 
SFINTtL,I)=DINT(I) 

105 CDNTINUE 
J=1 

106 CONTINUE 

ROT ( 1, 1 ) = ( OCOS ( PH 1 1 J > ) )**2 
RCT C 1 »2) = ( DS INI PHI ( J ) ) ) **2 
RCTI1»3)=.5#DSIN(2*PFIIJ)) 

ROT ( 1 , 4 )=0. ODO 
ROT 12, l J = RQT (1,2) 
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ROT {2 » 2 J -ROT £1,1) 

ROT t 2 r 3 )=— ROT (1,3) 

RCT { 2 ,4 )=ROT £ 1,4 ) 

RQTC3,l)=2*KGT (2,3) 

RCTC3,2)=2*RGTll,3) 

PGTt3,3)=0C0S(2*PHI( J) ) 

RCTf 3,4)=ROT(l,4) 

R0T£4,l) = R0TU,4) 

K0T£4,2)=RUT (1,4) 

ROT ( 4»3)=R0T (1,4) 

ROT £4,41=1. 00 0 
DC 122 K= 1 , 4 1 

IF t TtT A ( J 1 .GE. STETA1K1 .ANO. TETA(J) .LE. STETA t K+l) ) GOTO 123 

122 CONTINUE 
GGTG 500 

123 THETA £ 1 ,1)= SMATC K. 1 ) + £SHAT l K+l , 1 l-SMAT ( K , 11 )*(TETA£ Jl-STETAl K) ) / 
K STETA[K+11— STETA1KJ 1 

T HET A 1 2 » 21 = SMAT ( K , 2 )+ ISMAT l K+l, 2 1-SMAT ( K f 2) ) * I TETA £ J )-ST£TA< K) 1 / 
1 t STETA (K+l) -STETA £K) ) 

THETA( 3,31 = SMAT ( K, 3 1+ ( SMAT ( K+l ,3 l-SMAT 1 K ,3 ) > *£ TETA l J )-STETA£ K) ) / 
ItSTETAIK+l)— STETA£K) > 

IHETA£4,4)=THETA13,3) 

THETA £4, 31= SM AT ( K, 4 ) + £ SMAT { K + l ,4 l-SMAT ( K ,4 J )*£ TETA t J1-STETA£ K) )/ 
1(STETA(K + 1 )— STET AIK) 1 
THETA [ 3 i 4)=— THE T A ( 4, 3) 

THETN=. 5* ( T HET A{ 1, 1 1 +TFETA ( 2,2)1 
DO 160 1=1,4 * 

DC 161 J J= 1 ,4 

THETA( I ,JJ)=THETA( I ,JJ J/THETN 
161 CONTINUE 
160 CONTINUE 

DO 107 1=1,4 
SSF IN ( 1 1 = 0.000 
DC 147 K= 1 , 4 

S5FIN( I J=SSFIN(I)+ROT(I,K)*SFINT(l,K) 
sssr-im=ssFiNtn 

147 CONTINUE 

107 CONTINUE 

DC 170 1=1,4 

SF I NT t L * I )=SSSF III) 

170 CONTINUE 

DC 108 1=1,4 
SSFINi 1 1=0. ODD 
DO 148 K=l , 4 

SSFINI I)=SSFIN( l)+THETA( I ,K)*SF INT£ L»K) 

sssfh n=ssFiNm 

148 CONTINUE 

108 CONTINUE 
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DC 171 1=1,4 
SF It'll ( L, I )=SSSF II U 
171 CONTINUE 
NSC=NSC-1 
J = J + 1 

IF (NSC.GT.O) GOTO 106 
DC 109 1=1,4 

SFINTIL ,1 ) = SFINTIL, I )*IDEXP l-AS*GAMA) ) 

109 CONTINUE 

WR I TE (6 , 110 ) L 

1, RITE (6,1 11HSF INTI L, 1). 1=1,4) 

NTGT=L 
L— L+ 1 
GOTC 102 

501 DO 4002 N=1,NS 
DO 4003 M=1,NRA0 
SKAD=SRA(M)*S(N) 

NTQTT=0 
DG 5000 1=1,4 
5000 FINTl I ) =0. ODO 

DC 3939 L=1,NT0T 

IFIRADIL) -LE- 5RAD1 GO TO 4000 

GCTC 3999 

4000 NT OTT=N 10 TT+ 1 
DG 4001 1=1,4 

4C01 FINTt I)=FlNTtn+SFINTCL. I) 

3999 CONTINUE , 

WR I TE 1 6 ,4004 ) S<N),SRAIM) 

WRITE 16,98) NTOTT 

98 FORMAT l /T 10, ‘NUMBER CF PHOTCNS USED = »,I5) 

4004 FORMAT C//T10,*S VALU E= • , E5 . 2 , 5X , » S POT SIZE =*,F5„2) 
DC 112 1=1,4 

112 FINT(I)=FINTm/NTGTT 
NRITE(6,113) 

WR ITE l 6 , 1 141 (FINTII), 1=1.4) 

PCLR=(FINT( 1)-FINT!2))**2 
PGLR=P0LR+FINTI3)**2 
P0LR=PGLR+FINT<4)**2 
PGLR=SCRTI POLR) 

POLR=POLR/ I F INT ( 1 )+F INTI 2) ) 

WRITEC6*S9) POLR 
RFAC=FINTI2)/F1NT( 1) 

WRITE16,97) RFAC 

57 FCRMAT1/T1Q»*R— FACTCRIEY/EX) = *,F7.4,//) 

4003 CONTINUE 
4002 CCNTINUE 
500 STOP 
END 
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LISTINGS FOR POLYMIE AND DBMIE ROUTINES USED TO 
CALCULATE THE VOLUME SCATTERING FUNCTIONS 
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MAIN 

MAIN PROG POL YM I E( VECTOR ) 

THE FOLLOWING DOUBLE PRECISION INPUTS ARE REQUIRED: 

RF R=R£AL PART OF REFRACTIVE INDEX 
RF I=IMAGINARY PART OF REFRACTIVE INDEX 
RADU^UPP ERBOUND ON RADIUS (MICRONS ) 

WAVE=WAVELENGTH IN MICRONS 

All J=PARAMETER$ FOR DISTRIBUTION ONE 

B( I )=PARAMETE RS FOR DISTIBUTION TWO 

TH£TO( J I =VECT OR OF ANGLES FROM 0-90 (COMPLIMENTS ARE ALSO CALCl 
OTHER INPUTS ARE: 

JX = NUMBE R OF ANGLES FROM 0-90 

NRAD=NUMBER OF RADII BETWEEN O-RADU 

NPARA=NUMBER OF PARAMETERS IN DISTRIBUTION ONE 

NPARA2=N UMBER OF PARAMETERS IN DISTRIBUTION TWO 

TWO=LOGI CAL VARIABLE TO ENABLE THE USE OF TWO DISTRIBUTIONS 


TWO FUNCTION SUBPROGRAMS DI ST ( RAD* A ) AND DIST2(RA0,B) ARE REQUIRED 
IN ADDITION TO PDSMIE SUBROUTINE 


TWO DATA SETS (6 AND 8) ARE USED FOR OUTPUT; NORMALLY 6= PRINTER 
AND 8=TAPE 


10 FORM AT (3 D15. 5 ) 

11 F0RMAT12DI5.5 ) 

12 FORMAT (D15.5) 

13 FORMAT (D15.5»I5) 

14 FORMAT (215) 

15 FORMAT (L5 ) 

16 FORMAT (15) 

17 FORMAT (D15.5 ) ' 

20 FORMAT { 1H I I 

25 F0RMAT(//T10, 'ELEMENTS OF THE TRANSFORMATION MATRIX FOR A SPHERE 
1WITH SIZE PARAMETER = «,F15.5) 

30 F0RMAT(//T10, 'REFRACTIVE INDEX. REAL = • ,D15 .5 ,T60 , ' IMAG INARY » , D15 
1.5,//) 

35 F0RMAT(T3, 'ANGLE' ,T17,» SIGMA1* ,T31 , ■ SIGMA2 • ,T46 , » S IGMA3* ,T6I,'SIGM 
1A4ST76, * INTENSITY' ,T91, 'POLARIZATION'//) 

40 FO RMAT (F10.4,5E15»6,F15.4) 

45 FORMAT (//T10 , • EFFICIENCY FACTOR FOR EXT INCTION • ,E 15.6 ) 

50 FORMAT ( //T10 , * EFFICIENCY FACTOR FOR SCATTER ING* ,E 15 .6 ) 

55 FORMAT (//TiO » ' EFFICIENCY FACTOR FOR ABSORPTION' ,E15, 6) 

60 F0RMAT(//T10 t * ASYMMETRY FACTOR* ,E15.6//I 
70 FORMAT(//T10, • TOTAL TIME FOR THIS CASE IN SECONOS= *,F15.3//) 
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MAIN 

2 F0RMAK//T10, 'PROBABILITY FOR THIS SIZE PARAMETER = *,D15. 5,//) 

3 FORMAT(//T10, 'NORMALIZATION FACTOR FOR THIS SET OF SIZE PARA= ' , 
1D15.5,//) 

REALMS RFR,RFI»X»QEXT , OS CAT , QAB$*THETD( 100) ,PQEXT,PQSCAT ,PQABS 
80 F0RMAT(//T10, 'SCATTERING CROSS SECT ION* , E15.6 ) 

REALMS ELTRMX (4, 100, 2 ) , ALAM, CON, CTBRQS, A VCSTH, PELTMX t 4, 100 ,2 ) 
REAL *8 PAVCTH,THE( 100) ,PBSCAT 
REAL*4 AIN(100,2 ) , POLR { 100 , 2 ) 

RE AL *4 PAIN( 100,2) ,PPOLRI 100,2) 

REAL-^4 PAI(100,2),PP0L(100,2) 

REAL*8 PR0B2, PNQRM2 

REALMS PQEX,POSC A,PQAB,PBSCA, PAVCT »PELTM(4»IQ0»2) 

REAL*8 RADU, ORAD, WAVE, GAMMA, A< 20) ,PROB,B(20) 

LOGICAL WRN,TWO 
WRN=. FALSE. 

C0N=3. 1415926535 897932D+Q 
INTEGER NPARA.NPARA2 
90 READ 15,10) RFR ,RFI , WAVE 
READ (5,14) JX,NPARA 
READ (5,12) (THETD(I ) ,1=1 ,JX) 

READ (5,13) RADU,NRAD 
1 FORMAT (D15.5) 

00 5 1=1 , NPARA 
READ (5,1) A(I) 

5 CONTINUE 

PEAD(5, 15 ) TWO 
DO 95 1=1, JX 
95 THE ( I )=THETD( I ) 

IF (TWO) GO TO 61 
GO TO 62 

61 READ(5,16) NPARA2 
DO 62 I=1,NPARA2 
RE AD( 5 ,17) Bill 

62 CONTINUE 
PQEXT=O.ODO 
PQEX=0. ODO 
PQSCA=O.ODO 
PQAB=O.ODO 
PBSCA=O.ODO 
PAVCT=0. ODO 
PQSCAT=0.000 
PQABS=0. ODO 
PB SCAT=0 • ODO 
PAVCTH=0 .000 
DRAD=RAD U/NRAD 
DO 1000 J=1,JX 
DO 1000 K=l,2 
DO 999 1=1,4 
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PELTMXt I , J,K) =0.000 
PELTMtl, J,K)=0.0D0 
999 CONTINUE 

PAIN ( J,K )=0.000 
PA I ( J»K) =0.000 
PPOL (J,K )=0. ODO 
PPOLR(J»K)=O.ODO 
1000 CONTINUE 
RAD=0.0 
PNORM2=O.ODO 
IF (TWO) GO TO 91 
PN0RM2=1 .ODO 

91 CONTINUE 

PNORM=O.ODO 
TI ME 1=0. 0 
DO 3000 L-l, NRAD 
RAD=RAD+DRAD 
DO 100 J=1,JX 
100 THETD(J)=THE( J) 

X=2.3D0-C0N*RAD/WAVE 

PROB=DIST{RAD,A) 

IF (TWO) GO TO 63 
PR0B2=0. ODO 
GO TO 64 

63 PRQB2=DI ST2( RAD , B) 

64 CALL SETCLK 

CALL PDBMIE ( X ,RFR ,RFI , THETD, JX , QEXT , QSCAT , CTBRQS ,ELTRMX » WRN ) 
CALL READCL (TIME) 

IF (WRN) GO TO 1001 
PNORM=PNQRM+PROB 
PN0RMZ=PNQRM2+PRQB2 
TI ME1=TI ME1+T IME 
QABS=QEXT-QSCAT 
AVCSTH=CTBRQS/OSCAT 
DO 150 K=l*2 
DO 150 J=1,JX 

AIN( J,K> = ELTRMXd, J,KM-ELTRMX(2, J,K) 

POLR( J,KJ = { ELTRMX (2 » J , K)— ELT RMX ( 1 , J , K ) )/AIN(J,K) 

AIN ( J»K ) = ,5*AINl J,K) 

PAIN( J f K)=AINU,K)*PROB+PAIN(J,K) 

PAI (J,K)=AIN( J,K)*PR0B2+PAI ( J,K) 

PPOL ( J, K) =POLR ( J ,K)* PRGB2+PP0L ( J » K ) 

PPOLR ( J , K ) = PP0LR ( Jf K ) +PQLRI J »K )*PROB 
150 CONTINUE 

DO 2000 1=1,4 
DO 2000 J=1,JX 
DO 2000 K=l,2 

PELTMXt I , J,K)=PELTMXtI »J,K)+ELTRMX(I ,J,K)*PROB 
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MAIN 


original Pabp TC, 

op quaui^ 


PELTMI I, J, K)=PELTM< I, J,K)+ELTRMX( I , J,K)*PR0B2 
2000 CONTINUE 

WR ITEI6, 20 ) 

WRI TE (6 , 25) X 
WR IT E ( 6 f 30 ) RFR, RFI 
WRITE(6,35) 

WRITE(6,40> ( (THETDI J) , (ELTRMXII »J,1) ,1=1,4) , AI N ( J , 1 ) , PO LR ( J , 1 ) ) , 
1J= 1, JX) 

WRITE (8,40) ( {THETDI J ) »(ELTRMX( I , J , 1 ) , 1= 1, 4) , AIN (J , 1 ) , PO LRI J , I ) > , 
1J=1, JX) 

DO 200 J=1,JX 

THETDIJ) = 180. ODO—THETD ( J ) 

200 CONTINUE 
JMX=JX-l 
DO 210 J=1,JMX 
JJ=JX-J 

WRITE (6 ,40) (THETDI JO ) , I ELTRMX ( I , J J , 2 ) , 1= 1, 4 ) , AIN( J J. 2) ,POLR{ JJ,2) ) 
C WRITE 18,40)1 THETDI JJ ) , ( ELTRMXI I ,JJ, 2) ,1 = 1,4) , AI N( J J, 2 ) ,POLR ( JJ , 2 ) ) 

210 CONTINUE 

WRITE(6,45) QEXT 
WRITE(6,50) QSCAT 
WRITEI6.55) QA8S 
WRITE(6,60) AVCSTH 
WRITEI6.2) PROB 
WRITE (6, 2) PR0B2 
WRITEI6, 20) 

WRITE(6,70) TIME 
PQSCAT=PQSCAT+QSCAT*PROB 
PG)SCA=PQSCA+QSCAT*PR0B2 
PQEX=PQEX+QEXT*PR0B2 
PQAB=PQAB+QABS* PR082 
PBSCA=PBSCA+QSCAT*(RAD**2)*PR0B2 
PA VC T-P A VC T+AVCS TH* PR0B2 
PQEXT=PQEXT+QEXT*PROB 
PQAB$=PQABS+QABS*PRGB 
PBSCAT=PBSCAT+QSCAT*(RAD**2)*PR0B 
PAVCTH=PAVCTH+AVCSTH*PROB 
1001 WRN= .FALSE. 

3000 CONTINUE 

DO 4000 J=1,JX 
DO 4000 K=l»2 
DO 4001 1=1,4 

PELTMXII »J»K)=PELTMXII»J» K) /PNORM 
PELTMII, J,K)=PELTM( I , J , K ) /PN0RM2 
4001 CONTINUE 

PAIN! J,K) =PA INI J,K)/ PNORM 
PA I IJ,K)=PAI( J,K)/PN0RM2 
PPOL ( J,K) =PPOL I J , K )/ PN0RM2 


99 



o n 


' MAIN 


PPOLR(J,K)=PPOLR( J,K)/PNORM 

4J30 CONTINUE 
C END FILE 8 

PQSCAT=PQSCAT/PNORM 
PQEXT=PQEXT/PNORM 
PQABS=PQABS/PNGRM 
PBSCAT=P BSCAT^CON/PNORM 
PA VCTH=PA VCTH/ PNORM 
PQSCA=PQSCA/PNORM2 
PQEX=PQE X/PN0RM2 
PQAB=PQAB/PN0RM2 
PB SCA=PBSCA*CON/PNORM2 
PAVCT=PA VCT /PNORM 2 
DO 6000 J=l,JX 
6000 THETD(J)=THE{ J) 

WRITE (6,20) 

65 FORMAK//TIO, * ELEMENTS OF TRANFORMAT ION MATRIX FOR POL YDI SPERSI ON* 
1,//J 

WRITE (6, 650 
WRI TE ( 6, 30) RFR.RFI 
WRITE(6,35) 

WRITE (6, 40) ( (THETO(J) , (PELTMXt I » J ♦ 1 ) , 1=1*4) , PAI N ( J, 1) , P POLR ( J , 1 ) 
1), J=ItJX) 

WRITE (8,40) ( ( TH£TD( J ) , ( PELTMX (I » J »I ) ,1 = 1,4) , PAIN (J, U,PPOLR(J, I) 

1 ),J=1,JX> 

DO 5000 J=I,JX 
THETD( J) =180 .ODO-THETOIJ ) 

5000 CONTINUE 
JMX= JX— I 

DO 5001 J=l, JMX 
JJ=JX-J 

WRITE (6, 40) ( THETDt J J) ,( PELTMX (I , J J ,2) , I =1 ,4 > , PA IN ( J J, 2 ) ,PPOLR 
1 (0J,2) ) 

C WRITE (8,40) (THETD( JJ) , ( PELTMX ( I , J J , 2) , 1=1 ,4) ,PA IN l J J, 2 ) ,PPOLR 

C 1 ( J J , 2) ) 

5001 CONTINUE 

C END FILE 8 

WRITE(6»45) PQEXT 
WRITE (6, 50) PQSCAT 
WRITE(6, 55) PQA8S 
WRI TE (6 , 80) PBSCAT 
WRITE(6, 60) PAVCTH 
WRITE (6, 3) PNORM 
WRITE(6, 70) TIME1 
WR I TE ( 6 , 20) 

DO 5010 J=1,JX 
5010 THETDt J)=THE( J) 

IF (TWO) GOTO 5002 
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5002 

GO TO 5003 
WRITE (6,20) 

MAIN 


WR I TE (6 , 65 ) 
WRITEI6, 30) 

RFR, RFI 


WRITE{6,35) 
WR ITE (6, 40 } 

( {THETDI J) 


L) » J*1 » JX ) 


C 

WRITE (8, 40 ) 

I (THETD(J) 

C 

1) ,J=1,JX) 



00 5004 J*l, 

JX 


THETDI J) =18C 

i* ODO-THETD 

5 004 

CONTINUE 



JMX=JX-1 
DO 5005 J=l, 

JMX 


JJ*JX-J 
WRITE (6,40) 

(THETDI JJ ) 

C 

1 ( J J , 2 ) ) 
WRITE (8,40) 

{ THETDI JJ ) 

c 

1 ( J J , 2 ) ) 


5005 

CONTINUE 


C 

END FILE 8 



WRITE{6,45) 

PQEX 


WR I T E { 6 , 5 0 ) 

PQSCA 


WRITEC6,55) 

PQAB 


WRITE(6,80) 

P8SCA 


WR1TE{6,60) 

PAVCT 


WRITE(6,3) 

PN0RM2 


WRITE(6,70) 

TIME1 

5003 

WRITE{6,20) 

STOP 



END 



SJSfSSS 

, (PELTMU ,J,1) ,I*1,4),PAI(J,1) 
, (PELTM{I,J,1),I=1»4),PAI(J,1) 

( J) 

, (PELTM(I,JJ,2),I=1,4),PAI(JJ, 
,< PELTMU, JJ,2), 1 = 1,4) , PAH JJ, 
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2) ,PPOL 
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PDBMIE 


SUBROUTINE PDBMIE (X,RFR,RF t,THETD, JX ,QEXT,QSCAT,CTBRQS, ELTRMX, WRNt 
I) ‘ 

C RADIATION SCATTERED’ BY -A SPHERE. THIS SUBROUTINE CARRIES OUT ALL 

C SUBROUTINE FOR COMPUTING THE PARAMETERS OF THE ELECTROMAGNETIC 

C COMPUTATIONS IN SINGLE PRECISION ARITHMETIC. 

C THIS SUBROUTINE COMPUTES THE CAPITAL A FUNCTION BY MAKING USE OF 

C DOWNWARD RECORRENCE RELATIONSHIP. 

C X 0 SIZE PARAMETER OF THE SPHERE, C 2 * P I * RADIUS OF THE SPHERE > / 

C WAVELENGTH OF THE INCIDENT RADIATION!. 

C RFO REFRACTIVE INDEX OF THE MATERIAL OF THE SPHERE. COMPLEX 

C QUANTITY. .FORMO ( RFR - I * RFI J 

C THETD! J)0 ANGLE IN DEGREES BETWEEN THE DIRECTIONS OF THE INCIDENT 

C AND THE SCATTERED RADIATION. THETDCJJ IS OR = 90.0. 

C IF THETD(J) SHOULD HAPPEN TO 8E GREATER THAN 90.0, ENTER WITH 

C SUPPLEMENTARY VALUEO SEE COMMENTS BELOW ON ELTRMX.. 

C JXO TOTAL NUMBER OF THETD FOR WHICH THE COMPUTATION AREREQUIRDE. 

C JX SHOULD NOT EXCEED 200 UNLESS THE DIMENSIONS STATEMENTS 

C ARE APPROPRIATELY MODIFIED. 

C MAIN PROGRAM SHOULD ALSO HAVE REAL THETD { 200 ), ELTRMX ( 4,2 00 ,2 ) . 

C DEFINITIONS FOR THE FOLLOWING SYMBOLS CAN BE FOUND IN • LIGHT 

C SCATTERING BY SMALL PARTICLES, H. C. VAN DE HULST, JOHN WILEY + 

C SONS, INC., NEW YORK, L957 ». 

C QEXT82 EFFIEC IENCY FACTOR FOR EXTINCTION, VAN DE HULST, P.14 + 127 

C QSCAT82 EFFICIENCY FACTOR FOR SCATTERING, VAN DE HULST, P. 14 + 127. 

C CTBRQSO AVERAGE (COSINE THETA) * Q SCAT , VAN DE HULST, P. 128. 

C ELTRMXd , J,K) 0 ELEMENTS OF THE TRANSFORMATION MATRIX F »V AN DE HUL 

C ST, P.34,45 + 125. I = 10 ELEMENT M SUB 2.. I = 20ELEMENT M SUB 1.. 

C 1=30 ELEMENT S SUB 2L.. I = 40 ELEMENT D SUB 21... 

C ELTRMXII , J, 1 ) REPRESENTS THE ITH ELEMENT OF THE MATRIX FOR 

C THE ANGLE THETD( J ) . . ELTRMX ( I , J , 2 ) REPRESENTS THE ITH ELEMENT 

C OF THE MATRIX FOR THE ANGLE 180.0 - THETD(J) .. 

5 FORMAT ( 10X • THE VALUE OF THE SCATTERING ANGLE IS GREATER THAN 90.0 

$ DEGREES. IT IS * ,E15.4) 

FORMAT 1//10X' PLEASE READ THE COMMENTS*//) 

FORMAT ( // 10X* THE VALUE OF THE ARGUMENT JX IS GREATER THE 100*) 
FORMAT (// 10X * THE UPPER LIMIT FOR ACAP IS NOT ENOUGH. SUGGEST GET 
1DETAILED OUTPUT AND MODIFY SUBROUTINE*//) 

REAU-8 X, RX,RFR,RFI,QEXT,QSCAT,T(5),TA(4),TB( 2) , TC < 2) , TD { 2) , TE (2) , 

2 CTBRQS, ELTRMXl 4, 100,2 ), PI (3,100) ,TAU(3,100) , 

3 CSTHT( 100) ,SI2THT< 100), THETD) 100) 

COMPLEX* 16 RF,RRF,RRFX,WM1,FNA,FNB,TC1,TC2, WFN(2) , ACAP (8000) , 

2 FNAP, FNBP 
LOGICAL WRN 

9 F0RMAK//T10, 'WARNING, ACCURACY NOT ACHIEVED'//) 

TA ( 1 )Q REAL PART OF WFN(l).. TA(2)0 IMAGINARY PART OF WF N < 1 ) • . 

TA ( 3 ) 0 REAL PART OF WFN(2).. TA(4)0 IMAGINARY PART OF WFN(2).. 

TB1110 REAL PART DF FNA...TB(2)0 IMAGINARY PART OF FNA. . . 

TC { 1 ) 0 REAL PART OF FNB. . .TC (2 ) 0 IMAGINARY PART OF FNB... 
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TDI1J0 REAL PART OF FNAP . * TD(2) IMAGINARY PART OF FNAP... 

TE ( 1 ) 0 REAL PART OF FNBP... TE<2)0 IMAGINARY PART OF FNBP... 
FNAP *■ FNBP ARE THE PRECEDING VALUES OF FNA + FNB RESPECTIVELY. 
EQUIVALENCE IWFNU), TA 1 1 ) ) , IFNA, TB(i)>, I FNB, TC(D) 
EQUIVALENCE (FNAP, TD( 1) ) , (FNBP, TEH)) 

IF ( JX .LT. 101 ) GO TO 20 
WRITE (6, 7) 

WRITEI6, 6) 

STOP 1 

20 RF=DCMPLXIRFR ,-RFI ) 

RRF = 1.0D0/RF 
RX = l.ODO/X 
RRFX - RRF * RX 

T ( 1)=(X**2) s(RFR**2-»-RFI**2) 

r 1 1 ) =DSQRT (T(l) ) 

NMX1 = 1. 1000 * T( L) 

IF ( NMX1 .LT. 7999) GO TO 21 
WRITE (6, 8) 

STOP 2 
21 NMX2 - TCI) 

IF ( NMX1 .GT. 150) GO TO 22 
NMX 1 = 150 
NMX2 = 135 

22 AC AP (NMX 1 * 1 > = ( O.ODO, O.ODO ) 

00 23 N = 1, NMX1 
NN = NMX 1 - N + 1 

ACAP(NN) = INN+1) * RRFX - 1 .ODQ/ I (NN+i )*RRFX + ACAP I NN+1 ) ) 

23 CONTINUE 

DO 30 J = 1, JX 

IF ( THETD(J) .LT. O.ODO ) THETDIJ) = DABS! THETDI J ) ) 

IF I THETDIJ) .GT. O.ODO ) GO TO 24 
CSTHT(J) = l.ODO 
SI2THTIJ) = O.ODO 
GO TG 30 

24 IF ( THETDIJ) .GE. 90.0DO ) GO TO 25 

T(l) = I 3.1415926535897932D+0 * THETDI J ) )/180. DO 
CSTHTIJ) - DCOSI TI 1 ) ) 

SI2THTIJ) r l.ODO - CSTHTIJ)**,2 
GO TO 30 

25 IF I THETDIJ) .GT. 90.0D0 ) GO TO 28 
CSTHTIJ) = O.ODO 

SI2THTI J) = 1 .000 
GO TO 30 

28 WRITE 16, 5) THETDIJ) 

WRITEC6, 6) 

STOP 3 

30 CONTINUE 

DO 35 J = 1, JX 
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300 

35 


60 


65 


70 


PIU,J) = 0.000 
P I ( 2 f J ) = l.ODO 
T AlT( 1 » J ) = 0.000 
TAU (2, J) = CSTHT(J) 

CONTINUE 
T ( 1 ) = DCOS(X) 

T ( 2 ) = DSINIX) 

WM1=DCMPLX (T ( l) ,-T( 2) ) 

WFNt 1 ) -OCHPLX (T(2)»TU)) 

WFN ( 2 ) = RX » WFNU) - WM 1 
TCI = ACAP(l) * RRF + RX 
TC2 = ACAPU) RF + RX 

FNA = (TCI « TA ( 3 ) - TA(1)) / (TCi * WFN(2) - WFN(l>) 

FNB = (TC2 *• TA(3 ) - TA{ 1 ) ) / ( TC2 * WFNC2) - WFN(ll) 

FNAP = FNA 

FNBP = FNB 

T( 1) = 1. 5000 

TB(1) = T( I) * TBU ) 

TB(2) » T ( 1) * TB( 2 ) 

TC(l) - T ( 1 ) * TC(l) 


TC(2) » T ( I ) » TC( 2 ) 
DO 60 J =1,JX 
ELTRMXU, J,l> = TB( 1) 
ELTRMX (2 » J ,1 ) = TB ( 2 ) 
ELTRMX(3,J,1) = TC ( 1 ) 
ELTRMX(4,«J,1 ) = TC( 2 ) 
ELTRMXU, J, 2) = TB( 1) 
ELTRMXU, J, 2) = TB(2) 
ELTRMXI 3 » J , 2 ) = TC(l) 
ELTRMX14, J,2) = TC( 2 ) 
CONTINUE 


* P I ( 2 » J ) + TC ( 1 ) * TAUI2.J) 
+ PI(2»J ) +TC ( 2 ) * TAU ( 2, J ) 
4PI ( 2, J ) + TBU) * TAU( 2, J ) 

* P 1 ( 2, J ) + TB ( 2 ) * T AU( 2, J ) 

* PI(2,J) - TC (I ) * TAU(2, J) 

* PI ( 2, J ) - TC (2) * T AU ( 2, J ) 

* PI{ 2, J) — TB ( 1) * TAU(2,J) 

* PI { 2 » J ) - TB (2 ) * TAU( 2, J ) 


QEXT ■= 2.000 * ( TBU) + TCCU) 

QSCAT =(TB(L)^2 + TB(2)**2 + TC(1)*'U + TC ( 2 2) /3.75D0 


CTBRQS = 0.000 
N = 2 

T{I) = 2*N - 1 
T ( 2 ) = N - 1 
T ( 3 ) = 2 * N + 1 
0070 J=1 , JX 

PI( 3, J) = (T ( I)*=PI (2, J)4CSTHT( J)-N*PI ( l, J) )/T(2> 

TAU (3, J) = CSTHT(U)*MPI (3 » J ) — PI (1 , J ) )-T (1 )*S 1 2THT ( J >*P I ( 2, 0 )+TAU( 1,J 
1) 

CONTINUE 
WM1 = WFN(l) 

WFN(I) = WFN( 2 ) 

WFN(2) = TU) * RX * WFN(I) - WH1 
TCI = ACAP(N) * RRF + N * RX 
TC2 = ACAP(N) * RF + N -I- RX 
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FNA = (TCI ^TA(3) - TA(ll) / (TCI * WFN{2) - WFNdiJ 
FNB = ( TC2 * TA ( 3 ) - TA(U) / ( TC2 * WFN(2) - WFN(l)) 

T( 5 ) = N 

T(4I = T ( 1 ) / ( T(5 I * T( 2 ) ) 

T( 2) = ( T ( 2) 4 ( T (5 ) + 1.0D0))/T(5) 

CTBRQS = CTBRQS + T(2) <■' (TO(U * TB(1) + TD(2) * TB(2J + TEd) ' 
$TC(1J +■ TE { 2) * TC ( 2 ) ) + T(4) * (TD(1) * TE(1) + TD(2J T£<2)) 
QEXT = QEXT + T (3 ) * (TB(1) + TCdH 
T { 4) = TB(1)**2 + TB{ 2)'**'-2 + TC(l)+*-2 + TC(2)>r*2 
QSCAT = QSCAT + T(3) > T(4) 

T(2) = N (M + 1) 

T(l) = T ( 3 ) / T (2 ) 

K = (N / 2) * 2 
DO 83 J = I, JX 

ELTRMXd, J,1)=ELTRMX( 1, J , 1) +T{1 H ( TB ( 1 ) *PI ( 3 , J) +TC ( 1 )*TAU (3 » J) ) 
ELTRMXd, J,1>=ELTRMX( 2, J , 1 H-T { 1 )*(TB (2 )*P I ( 3, J H-TC ( 2 )*TAU( 3, J)) 
ELTRMX(3 , J ,1 ) =ELTRMX (3, J* I >+T (1 J -* (TC ( 1 )*P I ( 3, J )+TB ( 1 mAU( 3, J ) ) 
ELTRMX(4,J,1I=ELTRMX(4, J, 1M-T ( 1>* ( TC ( 2 )*PI ( 3 , J J + TB (2 >*=TAU (3 , J) ) 
IF(K.EQ.N) GO TO 75 

ELTRMX(1,J,2)=ELTRMX( 1, J , 2)*T{ 1 )* ( TB ( II *PI ( 3, J) -TC ( 1 ) dAU(3, J) ) 
ELTRMX(2,J,2>=ELTRMX(2» J,2J+T(1 l* (TB (2 ) *P I ( 3, J ) -TC ( 2 )*TAU (3, J I » 
ELTRMXd, J,2)=ELTRMX(3, J , 2H-T ( 1 ) * (TC (11 >rPI ( 3, J > -TB (1KTAU(3» J)> 

EL TRMX( 4 , J ,2 ) =ELTRMX (4 , J ,2 1 +T ( 1 )*( TC ( 2 ) *PI < 3, J J-TB ( 2 P-TAU ( 3, J ) 1 
G0T080 

75 ELTRMXd, J, 2) -ELTRMXd, J, 2) +T(ld(-TB( 1>*P 1(3, Jl+TCdWAU (3, J) ) 
ELTRMX (2, J,2 )=ELTRMX(2, J, 2)+T ( 1 ) * (-T8 ( 2) *P I ( 3,J l+TC ( 2) H TAU(3 , J) 1 
ELTRMXd, J, 21 =ELTRMX(3»J>2}+T{1)*(-TC{1)*PI ( 3 , J ) +T B( U^T AU ( 3, J )} 
ELTRMX(4»J,2) =£LTRMX ( 4, J» 2) +* ( 1 )d~TC { 2 )*P I ( 3, J 1 +TB( 2 )*T AU(3 , J) ) 
30 CONTINUE 

I F ( T (4 ) .LT. 1 .00— 14 ) GO TO 100 
N = N + 1 
DO 90 J = 1, JX 
PId, J) = PI (2, J) 

PI (2, JJ = PId, J> 

TAUd, J) •= TAU{ 2, J) 

TAU ( 2 , Ji = TAUd, J) 

90 CONTINUE 

FNAP = FNA 
FNBP - FNB 

IF (N .LE. NMX2 1 GO TO 65 
WRITE(6, 9) 

WRN = .TRUE. 

RETURN 

100 DO 120 J= 1 , JX 
00120K=1 , 2 
DO 1151 = 1, 4 
T ( I )=£LTRMX( I , J,K) 

115 CONTINUE 
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ELTRMX(2* J»K) = T<1>**2 * T{2)**2 
ELTRMX( 1 y J t K ) » T{3)**2 + T(4i**2 
ELTRMX(3, J,K) = T{1J*TC3> + TI2)*TC4) 
ELTRMXC4, J,K) = T(2)*TI3> - TC4»*Ttl) 
120 CONTINUE 

Til) = 2 • 000 * RX**2 
QE XT = QEXT * Til) 

QSCAT * QSCAT * Til) 

CTBRQS = 2.OD0 * CTBRQS * TCI) 

RETURN 

END 


106 



DIST 


FUNCTION DIST(RAD,AJ 
REAL*8 A(2O),RA0 
REAL*8 DIST'» B»C 
B=-A( 3) 

C=RAD**A(4> 

C=B*C 

DIST=A(1>*(RAD**A(2>>*DEXP<C J 

RETURN 

END 


ORIGINAL 

OF POOH QUAIJH1 
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DIST2 


FUNCTION DI ST2 ( R A0» B) 
REAL*8 B (20 ) »RAD 
REAL* 8 0IST2»A 
A=-(B(2»*1> 

DIST2-B(1I*B(2I*(RAD**A> 

RETURN 

END 


original 

*** < m ’6aSS 
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APPENDIX E 


PROGRAM LISTING FOR CURFIT ROUTINE USED TO FIT THE 
THEORETICAL SIZE DISTRIBUTIONS TO THE EMPERICAL DATA 
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MAIN 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE CURFIT 

MAKES A LEAST SQUARES FIT TO A NON-LINEAR FUNCTION 


DESCRIPTION OF PARAMETERS 

-ARRAY OF IND. VARIABLE DATA POINTS 
-ARRAY OF DEP. VARIABLE DATA POINTS 
-ARRAY OF STANDARD DEVIATIONS FOR Y 
-NUMBER OF DATA POINTS 


X 

Y 

SIGMAY 

NPTS 

NTERMS 

MODE 


DATA POINTS 


-NUMBER OF PARAMETERS 

-DETERMINES WEIGHTING FOR LEAST SQUARES FIT 
+ 1( INSTRUMENTAL) W(I )=1./SIGMAY( I)**2 
OINO WE IGHT ING )W( I 1=1 . 

-1 (STATISTICAL) W(I)=1./Y(I) 

A -ARRAY OF PARAMETERS 

DELTAA -ARRAY OF INCREMENTS FOR PARAMETERS 
FLAMDA -PROPORTION OF GRADIENT SEARCH INCLUDED 
YFIT -ARRAY OF CALCULATED VALUES OF Y 
CHISQR -REOUCED CHI SQUARE FOR FIT 


C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
C FUNCTN ( X, I * A) 

C EVALUATES THE FITTING FUNCTION FOR THE ITH TERM 

C SSP ROUTINE OSINV 
C INVERTS CURVATURE MATRIX 

C 

C COMMENTS 
C DATA FORMAT 
C NPTS»NTERMS,M0BE(3I5) 

C X(I),Y< n,(SIGMAY(I) ),(2(3>E12.6) 

DIMENSION X(100> ,Y( lOO), SIGMAY (100) ,AI20J , DELTAA (20) ,SlGMAAt20», 
IYPITUOO), YFITI (100) 

LOGICAL GRAD ,CUR,GRID 
21 FORMAT ( 3L5) 

R€AD(5»21) GRAD, CUR, GRID 
READ (5*1) NPTS, NTERMS, MODE 

1 FORMAT (31 5) 

IF (MCOE) 2,2,4 

2 READ (5,3) ( X( I ) , Y< I ) , 1=1 ,NPTS> 

3 FORMAT (2E12.6) 

GO TO 6 

4 READ( 5, 5 ) ( X( I ) , Y (I ) , SI GMAY ( I) , 1 = 1, NPTS ) 

5 FORMAT (3E12 *6) 

6 READ( 5, 7) ( A (J) , DELTAA( J } ,J = I , NTERMS ) 

7 FORMAT ( 2E12 »6 ) 

I SUM=0 
CHISQ1=1.0 

14 FLAMDA= . 0Q1 
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IF(CUR) GO TO 22 
I F < GRID J GO TO 23 

CALL GRADLS ( X»Y»SIGMAY,NPTS ,NTERMS,MGDE, A, DELTAA , 

1 YF I T »CHI SQR) 

GO TO 24 

22 CALL CUR FIT (X»Y,SIGMAY* NPTS, NTS RMS, MODE , A , DELTAA ,S IGMAA , FLAMDA, 
IYFI T,CHI SQR) 

GO TO 24 

23 CALL GRI DLS (X,Y ,SIGPAY»NPTS ,NTERMS »MODE»A,DELTAA, 
1SIGMAA,YFIT,CHISQR) 

GO TO 24 

24 PRINT 8, (A! J), J=1,NTERMS1 

8 FORMAT { • * » E 12« 6 ) 

PRINT 9, CHI SCR 

9 FORMAT C * ■ , *CHISQR=', 1X,E12.6,/) 

IF ( CHISQl— CHISQR ) 12,13,12 

12 CHISGi=CHISQR 
ISUM=ISUM+1 

IF (ISUM-10) 14,13,13 

13 DO 11 1=1 »NPTS 

11 YFITK I > = 1./YFIT(I ) 

PRINT 10 

10 FORMAT!* *,13X, 'INO.VAR.* , 12X, * DEP.VAR. * , HX , *1 NV. DEP .VAR .» , / » 
PRINT L5,CX{1), YFITtn.YFITHI} ,1=1, NPTSJ 
15 FORMAT! * « , 10X, E12. 6, 8X, E 12 .6, 8X ,E12 . 6 ) 

STOP 

END 
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CUR FIT' 

SUBROUTINE CUR FIT(X,Y,SIGMAY,NPT$fNTERMS , MODE, A 
1,DELTAA,$IGMAA,FLAMCA,YFIT ,CHSQR) 

DOUBLE PRECISION ARRAY 

DIMENSION X< IOO) , VI 100) , SIGH AY ( 100 > , A { 20 ) ,DELTAA (20) , SI GMAA { 20 > , 
1YFIT ( 100) ,W EIGHT! 100) , ALPHA! 20, 20) ,BETA( 20) » DERI VI 20 1 » ARRAY ( 20» 
120) ,B<20) 

11 NFREE=NPTS-NTERMS 
IF (NFREE) 13,13,20 
13 CHISQR=0. 

GO TO 110 

C EVALUATE WEIGHTS 

20 DO 30 1=1, NPTS 

21 IF { MCDE) 22,27,29 

22 IF ( Y( I ) ) 25,27,23 

23 WEIGHT! I ) = 1./Y( I ) 

GO TO 30 

25 WEIGHT! I )=l./<-Y( I) ) 

GC TC 30 
27 WE IGHTI I ) = 1 • 

GO TO 30 

29 WEIGHT! I )=1. /SIGMAY ( I >**2 

30 CONTINUE 

C EVALUATE ALPHA AND BETA MATRICES 

31 DO 34 J=1 , NTERMS 
BETA! J) =0. 

DO 34 K=1 , J 
34 ALPHA! J »K ) = 0. 

41 DO 50 1=1 , NPTS 

CALL FDE RI V ! X, I , A, DELTAA, NTERMS , DERI V ) 

DO 46 J=l* NTERMS 

BETA! J)=BETA( J)+WEIGHT ( I )* ! Y I I )-FUNCTN(X, I, A ) )*DERI V < J ) 

DO 46 K=1 » J 

46 ALPHA I J , K)=ALPH A ! J, K)+W EIGHT ( I l + DERIV ! J ) *DER IV(K ) 

50 CONTINUE 

51 DO 53 J=l, NTERMS 
DO 53 K=l» J 

53 ALP HA( K , J l=ALPHA ( J, K ) 

C EVALUATE CHISQR AT STARTING POINT 

61 DO 62 1=1, NPTS 

62 YFIT ( I)=FUNCTN{ X, I , A) 

63 CHI SQ1=FCHISQ!Y, SIGMA Y, NPTS, NFREE, MODE, YFIT) 

C INVERT CURVATURE MATRIX TO FIND NEW PARAMETERS 

71 DO 74 J=l, NTERMS 

72 00 73 K=1 , NTERMS 

73 ARRAY! J » K) = ALPH A ( J»K) /SORT ! ALPHA! J, J ) *ALPHA !K,K ) ) 

74 ARRAY! J,J)=1 -+FLAMDA 

80 CALL MAT I NV! ARRAY, NTERMS, 1) 

81 DO 84 J=l, NTERMS 
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CURFIT 


B( J)-A( J) 

DO 84 K=1,NTERMS 

84 JI+8ETAI Ki*ARRAY( J,K»/SQRT(ALPHAt J » J)*ALPHA{ K,K) J 

C IF CHI SQUARE INCREASED, I NCR EASE FLAMDA 

91 00 92 1=1 ,NPTS 

92 YFITC I J=FUNCTN( X ,1 , B) 

93 CHISQR=FCH1SQCY,SIGMAY,NPTS,NFREE,M0DE, YFIT> 

IF ICHISC1-CHISQR) 95,101,101 

95 FLAMDA=10.*FLAM0A ' 

GO TO 71' 

101 DO 103 J=l, NTERMS 
103 A(J|=B(J) 

FLAMDA=FLAHDA/10o 
110 RETURN 
END 

;KlGINAL PAGE'S? 

. »P POOR QUAIM 
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FDERIV 


SUBROUTINE FDERIVIX, I ,A,CELTAA,NTERMS t DERlV) 
DIMENSION X ( 100 ) »A( 20 J » DELT AA( 20 ) »DER IV ( 20) 
11 DO 18 J=l t NTERMS 
AJ=A { J) 

DELT A=DELTAA.( J) 

A(J)=AJ+DELTA 

YFIT=FUNCTN(X,I,A) 

A< J )=AJ-D£LTA ' 

DERIVC J J-( YF1T—FUNCTN(X» I»A) )/C2**DELTA) 

18 Al J) = AJ 
RETURN 
END 
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MATINV 


SUBROUTINE M ATI NV (ARRAY . NT E RMS , MCODE 1 
DOUBLE PRECISION ARRAY, B 
DIMENSION ARRAY (20,20) ,8(210) 

DO 1 1*1 , NTERMS 
DO 1 J=l,NTERMS 

CALL LQC ( I, J, I J ,NTERMS,NTER MS, MCODE ) 

1 B(IJ>*ARRAY(I,di 
EPS=1.0E-16 

CALL DSINV( B ,NT ERRS , EPS, IER > 

IF ( IERI 2,4,3 

2 PRINT 10 

10 FORMAT! * * , *N0 RESULT*, /I 
GO TO 4 

3 PRINT 11 

11 FORMAT! » •» ’WARNING* ,/) 

4 DO 5 1*1 , NTERMS 
DO 5 J*l, NTERMS 

CALL LOCU, J,IJ, NTERMS, NTERMS, MCODE) 

5 ARRAY(I,U)=B!IJ1 
RETURN 

END 


ORIGINAL PAGE I& 

° F poor quality 
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FCHISQ 


FUNCTION FCHISQ (Y, SIGMAY, NPTS*NFREE » MODE, YFI T> 
DIMENSION YI 100 J * SI GMAY 1 100 ) »YF IT ( 10.0 ) 

SUM-0 • 

DO 5 I=i,NPTS 
IFC MODE ) 1*2,3 

1 nn./mi 

GO TO 4 

2 W*l* 

GO TO 4 

3 W=1 ./ CS IGMAY II ) **2) 

4 SUM= { Y(I)-YFIT(IH*(Y{I I-YF IT! 1 1 i*W 

5 CONTINUE 
FCHISQ=SUM/NFREE 
RETURN 

END 
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SUBROUTINE DSINV 
PURPOSE 

INVERT A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX 
USAGE 

CALL DSINV(A»N,EPS, IER) 

DESCRIPTION OF PARAMETERS 

A - DOUBLE PRECISION UPPER TRIANGULAR PART OF GIVEN 

SYMMETRIC POSITIVE DEFINITE N BY N COEFFICIENT 
MATRIX. 

ON RETURN A CONTAINS THE RESULTANT UPPER 
TRIANGULAR MATRIX IN DOUBLE PRECISION. 

N - THE NUMBER OF ROWS (COLUMNS) IN GIVEN MATRIX. 

EPS - SINGLE PRECISION INPUT CONSTANT WHICH IS USED 

AS RELATIVE TOLERANCE FOR TEST ON LOSS OF 
SIGNIFICANCE. 

IER - RESULTING ERROR PARAMETER CODED AS FOLLOWS 
I ER=Q - NO ERROR 

I ER=— I - NO RESULT BECAUSE OF WRONG INPUT PARAME- 
TER N OR BECAUSE SOME RADICAND IS NGN- 
POSITIVE (MATRIX A IS NOT POSITIVE 
DEFINITE, POSSIBLY DUE TO LOSS OF SIGNI- 
FICANCE) 

I ER=K - WARNING WHICH INDICATES LOSS OF SIGNIFI- 
CANCE. THE RADICAND FORMED AT FACTORIZA- 
TION STEP K+l WAS STILL POSITIVE BUT NO 
LONGER GREATER THAN ABS(EPS A A(K+l,K+l) ) . 


REMARKS 

THE UPPER TRIANGULAR PART OF GIVEN MATRIX IS ASSUMED TO BE 
STORED COLUMNWISE IN N*(N+L)/2 SUCCESSIVE STORAGE LOCATIONS. 
IN THE SAME STORAGE LOCATIONS THE RESULTING UPPER TRIANGU- 
LAR MATRIX IS STORED COLUMNWISE TOO. 

THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL 
CALCULATED RADICANDS ARE POSITIVE. 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
DMFSD 

METHOD 

SOLUTION IS DONE USING FACTORIZATION BY SUBROUTINE DMFSD. 
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C 

SUBROUTINE DSINV(A,N, EPS, IER) 

C 

C 

DIMENSION A( 2 10 ) 

DOUBLE PRECISION A, DIN, WORK 

FACTORIZE GIVEN MATRIX BY MEANS OF SUBROUTINE DMF SO 
A = TRANSPOSE (T) * T 
CALL DMFSD(A,N,EPS, IER J 
IF(IER) 9,1,1 

INVERT UPPER TRIANGULAR MATRIX T 
PREPARE INVERSION-LOOP 

1 IPIV=N*(N+i)/2 
IND=I PIV 

INITIALIZE INVERSION-LOOP 
DO 6 1=1, N 
01 N=1.D0/A(IPIV> 

A( IPIV)=DIN 
MI N=N 
KEND= I— 1 
LANF=N-KEND 
IF(KEND) 5,5,2 

2 J= I NO 

INITIALIZE ROW-LOOP 
DO 4 K=1,KEND 
WORK=0. DO 
MIN=MIN-1 
LH OR=I PIV 
LVER= J 


START INNER LOOP 
DO 3 L=LANF,MIN 
LVER=LVER+1 
LHOR=LHOR+L 

3 WORK=WORK+A ( LVERl^A ( LHOR) 

C END OF INNER LOOP 

C 

A{ J)=— WORK*DIN 

4 J= J-MIN 

C END OF ROW-LOOP 

C 

5 IPI V=IPIV— MIN 

6 IND=IND— 1 

C END OF INVERSION-LOOP 
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CALCULATE INVERSES) BY MEANS OF INVERSECTJ 
INVERS E (A ) « INVERSE IT ) * TRANSPOSE ( I NVERSE( T) ) 
INITIALIZE MULTIPLICATION-LOOP 
00 8 1 = 1, N 
IPIV=IPI V+l 
J=IPIV 

INITIALIZE ROW-LOOP 
DO 8 K=I ,N 
W0RK = 0* DO 
LH0k = J 

START INNER LOOP 
DO 7 L=K,N 
LVER=LH0R+K-I 
wORKn=WORK+A(LHOR)*A(LVER) 

7 LHOR=LHOR+L 

END OF INNER LOOP 

A ( J ) =W0RK 

8 J=J+K 

END OF ROW- AND MULTIPLICATION-LOOP 

9 RETURN 
END 
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SUBROUTINE DMFSO 
PURPOSE 

FACTOR A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX 
USAGE 

CALL DMFSD(AjNjEPSflER) 

DESCRIPTION OF PARAMETERS 

A - DOUBLE PRECISION UPPER TRIANGULAR PART OF GIVEN 

SYMMETRIC POSITIVE DEFINITE N BY N COEFFICIENT 
MATRIX. 

ON RETURN A CONTAINS THE RESULTANT UPPER 
TRIANGULAR MATRIX IN DOUBLE PRECISION. 

N - THE NUMBER OF ROWS (COLUMNS} IN GIVEN MATRIX. 

EPS - SINGLE PRECISION INPUT CONSTANT WHICH IS USED 

AS RELATIVE TOLERANCE FOR TEST ON LOSS OF 
SIGNIFICANCE. 

IER - RESULTING ERROR PARAMETER CODED AS FOLLOWS 
IER=0 - NO ERROR 

IER=— 1 - NO RESULT BECAUSE OF WRONG INPUT PARAME- 
TER N OR BECAUSE SOME RADICAND IS NON- 
POSITIVE (MATRIX A IS NOT POSITIVE 
DEFINITE, POSSIBLY DUE TO LOSS OF SIGNI- 
FICANCE) 

IER=K - WARNING WHICH INDICATES LOSS OF SIGNIFI- 
CANCE. THE RADICAND FORMED AT FACTORIZA- 
TION STEP K+ I WAS STILL POSITIVE BUT NO 
LONGER GREATER THAN ABS( EPS*A <K+1,K+1 ) ) . 

REMARKS 

THE UPPER TRIANGULAR PART OF GIVEN MATRIX IS ASSUMED TO BE 
STORED COLUMNWISE IN N*lN+l)/2 SUCCESSIVE STORAGE LOCATIONS. 
IN THE SAME STORAGE LOCATIONS THE RESULTING UPPER TRIANGU- 
LAR MATRIX IS STORED COLUMNWISE TOO. 

THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL 
CALCULATED RAD1CANDS ARE POSITIVE. 

THE PRODUCT OF RETURNED DIAGONAL TERMS IS EQUAL TO THE 
SQUARE-ROOT OF THE DETERMINANT OF THE GIVEN MATRIX. 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 

METHOD 

SOLUTION IS DONE USING THE SQUARE-ROOT METHOD OF CHOLESKY. 
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OWGOJALPAOT^ 

OF POOR Q uALU 

THE GIVEN MATRIX IS REPRESENTED AS PRODUCT OF TWO TRIANGULAR 
MATRICES, WHERE THE LEFT HAND FACTOR IS THE TRANSPOSE OF 
THE RETURNED RIGHT HAND FACTOR. 


SUBROUTINE DM FSD( A, N , EPS, IER ) 


DIMENSION AI2I0J 

DOUBLE PRECISION DPIV,DSUM,A 

TEST ON WRONG INPUT PARAMETER N 
IF { N-l) 12,1,1 

1 IER=0 

INITIALIZE DIAGONAL-LOOP 
KP I V=0 
DO 11 K=1,N 
KPIV=KPI V+K 
IND=KPI V 
LEND=K— 1 

CALCULATE TOLERANCE 
TOL-ABS (EPS*SNGL(A(KPIV> ) ) 

START FACTORIZATION-LOOP OVER K-TH ROW 
DO 11 I=K,N 
DSUM=0. DO 
IF (LEND ) 2,4,2 

START INNER LOOP 

2 DO 3 L=1 , LENO_ 

LANF=KPI V— L 

LI ND=IND-L 

3 DSUM=DSUM+A(L ANF } ; fA (L IND) 

END OF INNER LOOP 

TRANSFORM ELEMENT A< INDJ 

4 DSUM=A( INDJ-DSUM 
IF(I-K) 10,5,10 

TEST FOR NEGATIVE PIVOT ELEMENT AND FOR LOSS OF SIGNIFICANCE 

5 IFISNGL ( DSUM )— TOL ) 6,6,9 

6 IF { DSUM ) 12,12,7 

7 IF C I ER I 8,6,9 

8 IE R=K— 1 
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DMFSD 


COMPUTE PIVOT ELEMENT 
DP IV=DSQRT CDSUM) 

At KPI V)=DPI V 
DPI V= I • DO/ DP IV 
GO TO II 

CALCULATE TERMS IN ROW 
A( IND) = DS UM*DPIV 
IND=I ND+I 

END OF DIAGONAL-LOOP 

RETURN 

IER=-I 

RETURN 

END 
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